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NUMERICAL STUDIES OF THREE-DIMENSIONAL 
BREAKDOWN IN TRAILING VORTEX WAKES 

By P. F. Evans and J. E. Hackett 
Lockheed-Georg ia Company 


SUMMARY 


Finite element, three-dimensional relaxation methods are used to calcu- 
late the development of vortex wakes behind aircraft for a considerable down- 
stream distance. The inclusion of a self-induction term in the solution, 
dependent upon local curvature and vortex core radius, permits calculation of 
finite lifetimes for systems for v;hich Infinite life would be predicted two- 
dimensional ly. The present report describes the associated computer program 
together with single-pair, twin-pair, and multiple-pair studies carried out 
using it. It is found, in single-pair studies, that there is a lower limit to 
the wavelengths at which the ”Crow"-type of instability can occur. Below 
this limit, self-induction effects cause the plane of the disturbance waves 
to rotate counter to the vortex direction. Self induction in two dimen- 
sionally generated twin spiral waves causes an increase in axial length 
which becomes more marked with decreasing initial wavelength. The time taken 
for vortex convergence toward the center plane is correspondingly increased. 
The limited parametric twin-pair study performed suggests that time~to- 
converge increases with increasing flap span. Limited studies of Boeing 7^7 
configurations show correct qualitative response to removal of the outer 
flap and to gear deployment, as compared with wind tunnel and flight test 
exper i ence. 
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1. INTRODUCTION 


During the past several years, strenuous efforts have been made to devise 
means for breaking up prematurely and permanently the trailing vortices behind 
large transport aircraft so as to reduce or eliminate their hazard to following 
aircraft. Many of the attempts aimed at modifying individual vortices have met 
with indifferent success or have been impractical or too costly to apply. 

One of the more reliable "natural” modes of vortex pair decay is vortex 
looping. This form of instability was analyzed first by Crow (ref. 1), using 
linear methods. Shortly thereafter, his results were extended numerically by 
Hackett and Theisen (ref. 2), who calculated the progression from Crow's 
sinuous shapes toward the familiar "looped" form sometimes observed in flight 
(ref. 3) or in water tank experiments (refs. 2 and 4). The calculations used 
finite-element, time-dependent techniques applied to single vortex pairs. A 
limited study was made of the effect of the initial amplitude and wavelength 
of a "Crow"-type perturbation. Because these initial quantities must be 
specified, there are difficulties In specifying time-to-loop. 

It is the objective of the present work to provide methods for deepening 
the previous single-pair studies and extending the methods to multiple-pair 
configurations including both wakes typical of current jumbo jets and experi- 
mental cases such as "sawtooth" loading suggested by Rossow (ref. 5) for 
promoting early wake disruption. Appropriate Fortran algorithms have been 
written and documented which achieve the above objectives. The program was 
delivered to NASA-Ames in April 1976 and implemented there on both the IBM 
360“67 and CDC -7600 computers. Program development was carried out at the 
Lockheed-Georg ia Company on a MAC-1 6 mi n i -computer , where CRT graphics are 
available and on a Univac 1106 terminal system for longer runs. The program 
and its documentation form Appendix A of the present report. 

in any study like the present one, it is very desirable to make compari- 
sons with experimental results. Though documented flow investigations in the 
near-field are quite numerous and despite the wealth of longer-time measure- 
ments at discrete downstream planes. It has proved impossible to find quanti- 
tative experiments which give adequate details of vortex looping flow. 
Experiments with towed models are the most suitable for this and although 
looping has been noted with a small model size relative to the test region 
cross section (ref. 2), the majority of scale models tested have had larger 
relative size and looping has not occurred. In some cases, wall constraints 
may have inhibited looping; in others the vortices have encountered a hori- 
zontal surface before looping had time to occur. Even with the looping flow 
properly generated, good measurements are difficult. 

In carrying out numerical calculations of the vortex wake, a choice is 
available between a space-fixed and an aircraft-fixed reference system. During 
the present work, It became apparent that while vortex looping is an unsteady 
space-fixed event, spiral multiple-vortex development is steady in an aircraft- 
fixed frame of reference. Though computing considerations made space-fixed 
analysis mandatory for a very long wake, it was found in practice that an 
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aircraft-fixed analysis was essential to provide proper initial conditions for 
multiple vortex arrays. An extended version of the program was therefore 
developed for this purpose, which Is compatible with the main, space-fixed 
algorithm. 

Sufficient experience has now been obtained with the program that the 
proper initialization procedures, step sizes, etc. are In hand and the signi- 
ficance of each input variable is well appreciated. These factors are docu- 
mented herein. However, since some of the quoted twinand multiple-pair 
examples lie part-way up the learning curve, the corresponding results should 
not be regarded as definitive but rather as a qualitative guide for future 
application of the program. 

Section 2 wi 1 1 outline the theoretical approach, velocity calculation 
methods, self-induction effects, stepping procedures, etc, used throughout the 
present work. Section 3 is devoted to single vortex pair events and further 
examines. the vortex looping phenomenon. In Section the distinction between 
looping and spiral center-plane convergence is discussed and the effects of 
twin pair variables are reviewed. The computational details and requirements 
for such calculations are also included. 

In Section 5 pilot calculations are described involving multiple vortex 
pair wakes for actual aircraft configurations. Conclusions are given in 
Section 6. 
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2. THEORETICAL APPROACH 


2.1 General Discussion 

When dealing with single-pair decay, there is little option but to impose 
an initial perturbation either arbitrarily or via a knowledge of a specific 
turbulence environment for aircraft flight. For multiple pai rs," however, this 
situation is relieved by the fact that in some cases vortex helices form 
naturally, starting in the near wake. In other cases, vortex strengths and 
positions may be such that the pattern disperses rapidly v/i thout helices form- 
ing. Figure 2.1 shows an example of each type, in end view, for twin vortex 
pairs, Bilanin er oL (ref. 6) have produced classification charts showing the 
relative vortex strengths and relative spans associated with each type of 
motion. 

in analyzing single-pair deformations, a space-fixed frame of reference 
IS perfectly adequate. However, it is evident from the above that an aircraft- 
fixed analysis is needed for multiple pairs, at least to start the calculation. 
Accordingly, a variant of the original program has been written, in an 
aircraft-fixed frame of reference, which includes the near-field effect of the 
wing. The resulting program solves the steady-state problem in three dimen- 
sions by a time relaxation of an initial two-dimensional, Trefftz plane 
solution. 

Once the initial, transient motions are complete there may remain a number 
of vortex pairs arranged in helices. it is not practical to continue to 
analyze the subsequent motion in aircraft-fixed axes, because of the prohibi- 
tive number of points involved, which must ail be perturbed starting from zero 
time. A change is therefore made Into space-fixed axes and a regular cyclic 
geometry is assumed. If this change is made at a distance corresponding to N 
miles aft of the aircraft and the subsequent time- to- looping is T seconds, the 
total distance to looping is given by (N + U„T) miles. 

It must be remarked that the above trans ient-pl us-cycl I c distinction is 
arbitrary, since divergence may require many cycles in some cases. However, 
the distinction is useful and practical except possibly for borderline, near- 
cycl ic cases. 


2.2 Velocity Calculations 


Given a constant-strength vortex, strength F, radius 6, which follows an 
arbitrary curved path in three dimensions from point A (x^, y/\, z/^) to point 
B (xg, yg, zg) , we seek the three velocity components, u, v, and w, at an 
arbitrary point N(x, y, z) on the vortex. l(x', y‘, z') is a typical element 
which influences N. 
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The calculation Is divided into three parts: the near, intermediate, and 

far fields. The near field concerns the effect of an element on itself, while 
the far field concerns elements sufficiently remote that their detailed struc- 
ture is not significant at point N. All other points are classified as 
i ntermediate. 

When studying the periodic vortex cases, it is convenient to use imaging 
techniques such that only one side of the aircraft wake and only one-half wave 
(one full wave for spiral perturbations) are actually perturbed and the 
periodic part of -the wake vortex system is represented by various reflections, 
in the case of aircraft-fixed analysis, there is no periodic geometry, but 
half-plane symmetry still applies. 


Near-field . — The self-induced velocity Is determined In a manner similar 
to that of Leonard (ref. 8 ). The self-induced velocity (Vs) is in a direction 
normal to the local plane of curva'ture and" Ts given by the relation 


Vs 




( 1 ) 


where: 


Rj. = local radius of curvature 

01, 02 = angles between the local point and each of the adjacent 
points based on the origin of the radius of curvature 

6 = core radius, derived from 


_ ii 

dt % dt 


(la) 


k 


Z = length of local segment 



The terms in the time-dependent expression for core radius, 6, yield 
respectively the viscous growth and vortex stretching effects. Since molecular 
viscosity has been used in the present calculations, viscous growth is very 
small. Equation (la) is evaluated in finite-difference form. 


Intermediate field . — We resolve the total strength F into orthogonal 
components F^, F 2 , and Fj. at x' , y', z' so that 


F^ + F^ + F^ = F^ 
X ^y ^ z ^ 


( 2 ) 


and 


F = 


dx'/ds' dy‘/ds' dz'/ds' 


(3) 


where 


dx''^ + dy'^ + dz''^ = ds 


I 2 _ Je I 2 


(i») 


Reference 7 gives: 
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So 
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B 

J 

A 




B 

f 1*7 

- — - (y -y') ds' 
J AiTr3 
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Writing Equation (5), at a point N, as a sum over LjOTAL finite elements: 


koTAL 


= I 


L=1 


(FAy) 


Ltotal 


(zn"Zl) - I 


L=1 


(TAz) / 

1 i" Wn -yl) 

HirrML^ 


■-TOTAL 

= I . , 3 

L=1 


(FAz) 


, , JTAXI, , 


> N 5^ L 


■-TOTAL 3 ■-TOTAL 


( 6 ) 


'"NL'- 
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where Ax, Ay, and Az are the components of a straight, finite element which 
represents a segment of the vortex AB. When N = L, Equation (l) is used, in 
implementing Equation (6) for N 5 ^ L, a "short element" concept is applied, in 
which the products, FAxl, FAyL, and FAz[_ are defined as the short element 
strength components. Dummy points {XDUM, YDUM, ZDUM) are situated midway 
between the perturbed points along each vortex line. Therefore: 

Axl = XDUMl+1 - XDUMl 

AX[_ = YDUMl+^ - LDUMl (7) 

Azl “ ZDUM[_+i - ZDUMl 


By using this short element method the effects of vortex stretching in the 
intermediate field can be considered, without changing the total number of 
points, by suitable modifications to FAx[_, FAy[_, and rAz[_. 

Far field . — Beyond one or two complete wavelengths, each side of a 
perturbed point the vortex system may be extended by horizontal, straight-line 
elements. These are added at both ends of periodic geometries (space fixed) 
and at the downstream end of aircraft-fixed geometries. 
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2.3 Time-History Calculation 


A one-step, second-order predictoi — corrector computation is used in which 
the desired element x-d isplacement, for example, is defined as the midpoint of 
the resultant of two Euler displacement vectors, u(t) dt and u(t + dt) dt. 
u(t) is the induced velocity field based on the positions of the vortex 
element system at the time t and u(t + dt) is the induced velocity based on the 
positions of the system at time (t + dt). The general relation for displacement 
using predictor-corrector stepping is: 

x'(t + dt)=y[x(t)+x{t + dt)+u(t + dt) *dt] (8) 


where 

x(t + dt) = x{t) + u(t) • dt (Euler value at (t + dt)) 
x(t) = value at t 

x'(t-i-dt) = predictor-corrector value at t + dt 
The other two components are treated similarly. 

in two-dimensional calculations the accuracy of the t ime- i ntegrati-on may 
be checked fairly readily by determining the vortex center of gravity, moments, 
and the Ki rkof f-Routh path function (see ref: 5)* Such checks are generally 
not available for three-dimensional motions. In principle, checks should be 
made that conservation laws are obeyed, but the necessary volume integrations 
are more difficult and more numerically hazardous than the original calcula- 
tion. The only check found feasible so far is to employ a reversed-march i ng 
technique to find the degree of accuracy to which the system returns to its 
original state. At the end of a series of time steps, the sign of the time 
increment is reversed and the program marches back to the initial time. The 
deviation for each vortex element is then defined as; 


21 ^ 


d| = [(x{f - Xj^)2 + (yj^ - yj^)2 + 


(9) 


where f denotes position after backward marching, and o denotes initial 
position. The arithmetic 


N 


I ■'i) 


and geometric 


i=l 


(- (dj 62 df^[) ^ ) 


( 10 ) 


Some typical results are illustrated in the 


means are then determined, 
follow! ng, table: 
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TABLE 1. - ACCURACY OF TIME-STEPPING METHODS 
(SINGLE PAIR, 51 STEPS) 


Calculation Method 


Arith, Mean 


Geom. Mean 


Predictor-Corrector .00313 

Euler .801 


.00271 

.696 


The predictor-correlator method, which takes little more than twice the running 
time of the simple Euler stepping method, provides a 250-fold accuracy gain in 
this particular example. In practice, a part of such benefits is used to 
reduce computer run time. 


2. A Time-to-Loop 


The formation of long loops along the track of an originally parallel 
trailing vortex pair is a well-known occurrence. The development from the 
initial sinuous perturbation (ref. 1) towards the near-ring shapes may be 
calculated (ref. 2) by finite element methods. The next stage, the formation 
of loopSjis not well understood. In real flows, the point of nearest approach, 
across the plane of symmetry, starts to convect downwards; but almost immedi- 
ately there is a "click" and the vortices break, change partners, and then link 
across the center-plane to form loops. This event is a very fundamental one 
since it also typifies the cascading process by which large eddies break down 
to smaller ones in turbulent flows. Nevertheless, there is Kittle detailed 
knowledge of the mechanism, and no criterion exists which defines when the 
"click" occurs. 

Fortunately, it is sufficiently accurate for present purposes to define 
the start of looping as the time at which the point of closest approach is 
convected downward at ten times the mean downdrift velocity for the whole wake. 
This defines the start of looping sufficiently sharply- for use in long-time 
solutions. (See Figure 2.2.) 
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3. SINGLE-PAIR LOOPING 


3.1 Significant Parameters 


Reference 2 gives examples which extend numerically the linearized analysis 
by Crow (ref. 1). As indicated above, these space-fixed calculations have been 
improved, both numerically and by the addition of the element self-induction 
term of Equation (I). As a result, the tendency to a fine-grain waviness in- 
stability on each individual loop has been reduced very considerably. The use 
of the higher preci s ion- CDC- 76 OO computer for many of the runs reported here 
.gave additional numerical improvements and permitted longer time histories to be 
calculated. Figures 3.1(a), (b) , and (c) show, respectively, plan views of 
the initial condition, a later geometry using results from a reference 2 type of 
calculation and the corresponding result from the present program. 

For a given vortex spacing (b) , the time to loop (tj^) , as defined in Section 
2.k, is a function of the vortex strength (F) ; the amplitudes (a^^, ay), phase (B) , 
and wavelength (Xb) of the Initial perturbation and the vortex core radius (6). 

The relationships take the nondimensional form: 



3.2 Effect of Perturbation Amplitude 


In Section 3, -the initial perturbations and ay usually will be equal 
and in phase, giving the classical 1}5°- 1 ncl i ned plane wave shown in Figure 
3 . 1 (a). in Figure 3.2, it is evident that there is an approximately inverse 
relationship betv/een time-to-loop and initial amplitude. The constant of pro- 
portionality is, of course, dependent upon vortex core radius and wavelength. 

As perturbation amplitude increases, the drift velocity of the point-of-close- 
approach will increase. Eventually, this will reach ten times the mean 
(straight vortex) drift velocity as used to define the onset of looping. This 
produces the intercept with the horizontal axis. In potential flow, this would 
correspond approximately to the condition when the distance of closest approach 
was 10^ of the mean vortex spacing. 

Limited checks have also been made of the effect of using spiral and 
planar initial conditions. Labelled points in Figure 3-2 show that the use of 
a spiral Initial condition increases the time-to-loop by nearly k0% . This is 
important because multiple vortex pairs develop spiral modes. Clearly, the 
"Crow" mode is more efficient at driving the system towards looping. It is 
Interesting to note that, if the x- and y-perturbat ions are applied separately, 
the x-perturbation on its own is very inefficient in this regard. However, 
applying solely y-perturbation increases time-to-loop by only 20% over the 
combined perturbation. 
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These results suggest that the present methods could be used to determine 
optimum configurations for deliberately induced initial perturbations in real 
f 1 ows . 


3.3 Effect of Vortex Core Radius 

Figure 3-3 shows the results of three calculations made for the same total 
time. The first uses current numerics, but omits the self-induction effect on 
each element. Figures 3-3Cb) and (c) show corresponding results for medium and 
large-sized cores. As might be expected, increase in core size has a damping 
effect on the overall motion. It is also interesting to note that the omission 
of curvature effects produces a result intermediate to the other two; i.e. an 
effective core radius has been implied. This is a somewhat similar situation 
to the use of a cut-off radius by Crow (ref. l). In the present case, the 
effective radius is determined (If self-terms are neglected) by the number of 
elements per wave. 

Figure 3-^ returns to the full equations and shows a gradual increase in 
time-to-loop with vortex core radius. It should be noted that the right hand 
part of this curve becomes less valid because the small radius assumptions in- 
herent in Equation (1) are violated. 

A number of Investigators have suggested means for diffusing individual 
vortices to reduce their hazard to following aircraft. The present result 
suggests that such diffusion may actually increase the time to ultimate decay, 
by delaying looping. It is therefore important to be certain that such 
attempts to accelerate diffusion have a sufficient effect. 


’ 3 .^ Effect of Perturbation Wavelength 


Though time-to-loop decreased with decrease in wavelength, as previously 
(ref. 2) , an unexpected upturn in the curve occurred at short wavelengths in 
the present calculations (see Figure 3-5)- 

To better understand the mechanisms, let us first examine in more detail 
the progression towards "conventional" looping. Figure 3-6 shows, at zero 
time, an end view of the right-hand vortex, incl ined outward at ^5 degrees. 

The individual points are calculation points on a half-wave. As time pro- 
gresses the whole array drifts downward, largely maintaining its 45~degree 
inclination. At the two longest times shown, the lower points are moving 
downward and inward at an increasing rate: a nascent looping stage has been 

reached, 

I^igure 3.7 has essentially the same conditions as previously, but a 
shorter wavelength. Initially, the innermost vortices move inward as before 
but at a faster rate, as would be expected from previous experience concerning 
wavelength effects. However, the 45“degree inclination is not retained as 
before. A clockwise rotation, opposing the direction of the vortex, develops. 
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In consequence, the di stance-of-closest approach starts to increase markedly 
and looping is forestalled. 


The suggested mechanism for the above motion is illustrated in Figure 
3 . 8 (a). At any wavelength, there is a tendency for the crests of each sine wave 
to self-convect like part of a vortex ring. The corresponding motion is out- 
ward and downward at the top of the wave, inward and upward at the bottom, and 
zero at the inflection points. This motion increases as the wavelength is 
reduced. Counterbalancing th] s effect is the velocity gradient due to the 
opposite vortex, which Is approximately constant. Thus, the outward rotation 
due to self-induction, which is dependent upon the ratio of amplitude to wave- 
length, is offset against the velocity gradient effect, which depends only upon 
amplitude. An illustrative velocity vector diagram is shown in Figure 3.8(b). 

Boundaries have been found, in an ampl i tude ~ wavel ength domain, which 
define the lower limit for looping. Figure 3-9 shows the results of a number 
of numerical calculations, for two core radii. As might be expected from cored - 
ring convection properties, the larger vortex core radii cases are less likely 
to self-rotate. 

To each point on the boundaries shown in Figure 3.9> there corresponds a 
(different) time-to-loop. These are minima, as was noted in Figure 3*5- 
Figure 3.10 illustrates the nature of time-to-loop contours in the same 
coordinates as the previous figure. 
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4. TWIN PAIR CONVERGENCE 


4,1 The Generation and Disruption of Spiral Waves 


It was noted previously that single spiral waves were noticeably more 
stable than plane waves of the "Crow” type. In this section we shall examine the 
properties of spiral waves with emphasis on the mechanisms and flow geometries 
which lead to vortex destruction. Though the same calculation procedures and 
divergence criteria may be used for each, there is a very important distinction 
between the instability of a spiral wave and that of a plane wave of the ”Crow”- 
type. Figure 4.1 illustrates this distinction from the viewpoint of an observer 
fixed in space. When a vortex trail with plane waves reaches a sufficient age, 
this observer sees "necking” start to occur, followed by the aforementioned 
"click" when looping occurs. Thereafter, the looping point remains fixed in 
space, IP contrast, the spiral wave is fed continuously towards the center 
plane. Prior to the critical time, the fixed-observer sees a convergence event 
moving towards him in the direction of aircraft motion. The "time-to-converge" 
corresponds to its passage directly in front of him. In contrast to the 
previous case, the event is steady as seen from an aircraft frame of reference. 
These considerations in no way preclude space-fixed events due to atmospheric 
disturbance or other modes of instability. However, to under! fne the distinc- 
tion, the distraction of spiral waves by meeting at the center plane will be 
referred to as "convergence" and "looping" will be reserved for the space-fixed 
event. ' 

Definition of variables . The spiral wavelength for a twin-pair trail is 
determined by aircraft forward speed, as well as vortex details and it must be 
included in the normalization scheme. Thus, if primed quantities are dimen- 
sionless and the kinematic variables are as shown in Figure 4.2, we define 


r-' = — 
'F U 


b/2 


T ' s * 

U b/2 

CO 


with a total 

For convenience, we write 


<|) = Tp* + Tj* 

- £l_ ’^t' 

^ Tp Tp' 


The following 


relationships are noted: 

p_‘ = — ^ — pi = — 

F 1 +Y ^T 1 +Y 

np + YHj = (1 +y)ticg' 


( 12 ) 

(T3) 


(14) 


(15) 

( 16 ) 
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( 17 ) 


We normalize time t (secs), after aircraft passage, via 

t U 

CO 

" b/2 

For multiple vortex arrays and because of the possibility of c.g. variations 
under viscous, three-dimensional conditions, we measure wavelength of spiral 
motions, etc., in aircraft spans, as Ab, it may be shown that the time per 
revolution of rp and Fy about their c.g., in two dimensions, is given by 

x/rev. = (r,.^ - ' (ig) 


or in aircraft variables 


o (hc + Yr)^) 

x/rev. ^ ^ ^ ( 19 ) 

and we note also that the wavelength, 

A = “ t/rev. • (20) 

Conditions for convergence . Figures 1t.3 and 4.4 show the effects, on the 
calculated development towards convergence, of varying the relative span posi- 
tions and strengths of the vortices. The vortex wavelength, after three- 
dimensional relaxation, is a function of these variables and both end view 
and t-ime-to- loop follow from this. Vortex radius also becomes a significant 
variable. 


With a wing tip vortex somewhat weaker than the flap vortex, Figure 4.4 
indicates that small flap spans lead to early convergence of the tip vortex. 

As the flap span increases, the spiral amplitude is reduced, and the trajec- 
tories approach the center plane less closely. Consequently, convergence Is 
delayed. A simple geometric analysis shows that when the flap vortex Is the 
stronger, the tip vortex will intersect the center plane after l80-degrees or 
less of travel provided that 

np < , (y £ 1). (21) 

Equation (21) Is a particular case of a geometric expression for minimum 
centerplane clearance by the tip vortex, i.e, 

^MiN = (Tip + y) - 1 , (y < 1) . (22) 

Figure 4.4 illustrates the consequences of Equation (22). With flap-span 
constant, a case with y less than one [Figure 4.4(a)] clearly converges much more 
quickly than with equal strength [Figure 4.4(b)]. We note here that since the 
two spirals are of equal amplitude, they each have an equal chance of converging. 
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Once flap vortex (now weaker) forms the outermost spiral and conse- 
quently clearly cannot be less than np. This explains the range of 

validity of Equation (22). 


The velocity condition for convergence . Equations (21) and (22), being 
based solely upon geometric considerations, do not reflect the velocity condi- 
tion used to define convergence in the present calculations. The effect is to 
move the critical point outward from the center plane by an amount 


y(t1f +y) 

N(Y + 1)2 


(23) 


which is added to Equation (22) to give 


n 


1 




MIN 1 +Y 


■} - 1 


Y < 1 


(24) 


S i nee y ^ 1 and an 
evidently smal 1 . 


N-value of 10.0 Is used, the effect of the added term Is 
A corresponding analysis modifies Equation (21) to the form 


’’F - (2(v^!)N^ y) ' ' 

This reduces to Equation (2l) for N->“ or Y“^0. Again, the differences are 
small within the range of validity. The differences between results derived 
from (24) and (25) and those from the simpler expressions are not sufficient 
to justify the added complexity in most applications. 

The spiral strain index, a . Equations (21) and (25) describe the condi- 
tions under which the flap vortex drives the tip vortex Into, or near to the 
the center plane after half a revolution or less of travel. The tip vortex, 
being the weaker, forms the outermost spiral, which converges. However, it 
Is evident from Figures 4.3 and 4.4 that, for cases with several revolutions 
prior to convergence, it is also possible for the Inner, more intense, spiral 
to converge first. Figure 4.5 shovis a further example of each type. 

Figure 4.6 Is a parametric map with flap span as ordinate and relative 
vortex strength as abscissa. Each plotted point depicts one computer run in 
a parametric study which will be described in Section 4.2. Triangles and 
circles represent, respectively, flap and tip vortex convergence as the first 
event. As might be expected, tip vortex convergence o'ceurs in the vicinity 
of the lower-left boundary [Eqns. (21) and (25)], but then flap vortex con- 
vergence becomes more likely as np is increased. There is also a change on 
passing through y= 1.0 (i.e. flap and tip strengths equal). To understand 
the structure of these domains, we shall examine the distortion mechanism 
which leads ultimately to convergence. 


On passing through its Innermost position, each vortex receives an extra 
downward impulse Induced by the opposite-hand vortex at a distance 2riu|f^ away. 
This causes an outward-spiraling tendency which, continued long enough, leads 
to convergence at the center plane. We may quantify this tendency in terms 



of the nondi mens Iona 1 velocity increment causing the effect. Thus, we define 
sp i ral -strai n indices: 


and 


r ' 


^MINF 


= 


ill 

’^MINT 


( 26 ) 


We now hypothesize that the vortex which has the larger a value vyill converge 
first. Thus, for the tip vortex to converge first requires that 


0T > ap 


Using Equation (22) and noting that yields 


rp‘y(l +y) Tp' 

^’^'^CRIT ~ ’^fCRiT 


(27) 


for the tip vortex to converge first. 

A borderline case may be found by using an equality in Equation (27) 
which readily yields 


•"CRIT y + 2 


( 28 ) 


This line partitions the data in Figure h~5 quite successfully though not 
perfectly because of the simplifying assumptions. However, resolution of the 
flap- or tip-convergent regions requires some care. 

Provided that ^Fq^ij ^f’ ^re both positive, Equation (27) may be 
wr I tten 


li 

D 


•^CRIT 

^’^'^CRIT 


y(l +y) 

- (1 -y) 


> 


1 


We may clear this further, by multiplying by D, only if 


i .e. if y > 1 - 2np^j^j.p. 

This defines a region In Figure 4.5 above the line marked "Eqn. (2l).‘‘ 

Thus , 


15 





giving 


(Y^+y-Z) > (y-1), Y> (1 -2^Pj,j^l^) (29) 

Since both sides change sign at y = l, the cases y>l must be considered 
separately. This yields, eventually 

'^f'CRlT ^ 7^ ’ 

and 

Dp - < — 3 

FcRIT y + 2 ’ 

for the tip vortex to converge first. This agrees with most of the test cases 
in Figure 4,6. More elaborate treatments are obviously possible, using Equa- 
tion (25) j for example, or taking spiral diameter into account. However, the 
simplicity of Equation (30) is attractive. 


4.2 Twin-Pair, Parametric Study 

By suitable choices of flap span and loading. It is possible to reduce 
the vortex hazard for any specific aircraft at a given weight and speed. 

Recent studies (ref. 11) have demonstrated the hazard reduction benefits of 
reduced flap span, for example. In order to provide a basis for quantifying 
the trades between flap span and setting angle, convergence times have been 
calculated over a wide range of rip and Fy/Tp. Boundaries for rapid conver- 
gence have been published previously (ref. 6), but no attempt has been made 
to calculate the time taken for cases within such limits. This Is only 
possible if three-dimensional, self-induction effects are included. 

Figure 4.7a shows contours of the dimensionless time-to-converge [see 
Equation (17") ]• for computational consistency a constant waverlength 
X = 5*375> was employed for the^ whole figure, the implicit assumption 
being that X is a weak parameter. However, subsequent studies (see Section 
4 . 3 ) have shown this to be untrue. Consequently, the loading parameter 
(C[_/PR) cannot be merged with nond i mens Iona 1 time t, as was originally hoped, 
it follows that a map like Figure 4.7a should be prepared for each of a series 
of constant values of Ci^/PR. In the present case the relevant Cl/AR value for 
any specific point may be found from Figure 4.7b using the fact that X =5-375 
in the previous figure. Figure 4.8 presents the results of a limited study 
(for Cl/AR = 0 . 25 ) using proper wavelengths for each run. An increase in tlme- 
-to-converge with increasing flap span is indicated, in view of the signifi- 
cance of other parameters (notably vortex radius), there is a very real 
question whether it is preferable to treat even twin pair cases on an 
individual, rather than a parametric basis, 
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Y>1 


(30) 





When convergent vortices reach the center plane, they merge and dissipate. 
However, it is possible and even probable, in still air, that the surviving 
vortex will still represent a hazard. Figure ^.3 shows that, for short con- 
vergence times near the lower hatched boundary, the surviving vortex is un- 
desirably strong. At any given flap span (>.333b) it is evident that the 
surviving vortex strength is maximized if flap- and tip-vortex strengths are 
equal. It Is possible to adopt a systems approach to define a minimum hazard 
condition and then determine optimum flap spans and vortex strengths for a 
given flight condition with both the surviving vortex strength (Figure 4.3) 
and time-to-converge (Figure 4.7a) taken into account. However, Figure 4.7a 
should be revised to include values of spiral wavelength consistent with the 
relevant C|_/AR condition before such studies are undertaken. 


4.3 Three-Dimensional Effects 

As indicated above, changes in spiral wavelengths were found to have a 
significant effect on time-to-converge. Figure 4.10, which was calculated on 
the same space-fixed basis as the previous study, shows an increase in time- 
to-converge as wavelength is increased. When two-dimensional calculations 
were attempted, to provide a limiting case, convergence was not obtained even 
well beyond t values corresponding to the six-mile limit generally applied. 

This indicates agreement between the present numerical result and Bilan in's 
(ref. 6) analytical result which predicts infinite time for two-dimensional 
cases outside the lower boundary. 

The detailed reasons for wavelength dependence are not understood. How- 
ever, it is noted that the self-convection effect on an element is similar In 
nature in a spiral to that found previously for plane waves; there is a 
tendency to slow down the gross rotation rate. Another distinction concerns 
“axial pumping." A tightly-wound spiral vortex has transverse vorticity com- 
ponents, which correspond to “ring" vorticity. In the core of the spiral, 
these vorticity components Induce a velocity away from the generating wing and 
a proportion Is felt by the spiral Itself, tending to lengthen it. It follows 
that any vortex spiral geometry generated on the basis of quasi-two-dimensional 
assumptions will have too short a wavelength. Figure 4.11 quantifies this 
effect. Vortex spirals were first calculated two-d imens Iona 1 ly , then relaxed 
three-dimensional ly In a wing-fixed frame of reference. In view of the X- 
dependence of time-to-converge (Figure 4.10), it is clearly desirable to 
determine a 3~D relaxed wavelength, rather than one derived simply from two- 
dimensional considerations, before embarking on detailed or definitive 
stud i es . 


4.4 Notes on Computing 


initial i zation . We saw in the preceding section that a spiral wavelength 
assigned on the basis of two-dimensional considerations [Eqns. (13) and (20)] 
will not be in equilibrium axially. It has also been necessary to assign unit 
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value to the tip vortex span, in the Interest of reducing the number of vari- 
ables. To circumvent such difficulties, it is now clear that any space-fixed 
calculation involving more than one vortex pair should at least start with a 
corresponding aircraft-fixed calculation which determines vortex spans and 
wavelengths, which are in self-equilibrium and which are consistent with the 
span loading. The two-dimensions-plus-time type of calculation provides a 
suitable starting geometry for three-dimensional relaxation of an appropriate 
(i.e. more than one-full wavelength) length of wake. 

Since it is unfeasible to relax the equivalent of several miles of wake 
simultaneously, in aircraft-fixed axes, a transition to a single-wave, space- 
fixed calculation is made. A suitable whole wavelength is selected* which is 
downstream of the direct influence of the wing; this is assigned time value 
T=0 in the space-fixed calculation. Figure 4.12 shows how subsequent time 
steps march down the wake. There are, thus, three computational components 
to the distance of any particular point from the wing. These are: X[ , the 

distance from the wing to the leading edge of the frame t = 0; t, which defines 
the distance travelled by the leading edge of the frame; and finally xq, the 
axial distance within the frame t, measured from its leading edge. 


The trajectory of a particular point . Since the point which converges 
is of predominant interest In the present study, we shall use this as an 
example. At the start of the calculation, it cannot be known which point will 
converge. Consequently, all points on every frame must be output. Once the 
convergent point has been identif.ied, however, its complete trajectory may be 
determined, frame-by-frame, back to the point Xq on the initial frame. The 
examples in Figure 4.5 were plotted in this way, Tq, which corresponds to the 
position on the original frame, is given by 


2N - 1 
NTOT 


2X 


(31) 


where N is the vortex element number and NTOT is the number of elements per 
wave. The total time-to-converge is thus given by 

cc = 1 + ■To + Tc 

r 

where T| equals X|/U^ and t^, as before, is the time step (space fixed) 
which convergence occurred. Any other specific point may be treated in 
same way. 

Conjugate runs . In twin-pair parametric studies, advantage may be taken 
of the fact that, after the first half wavelength, the flap and tip vortices 
Interchange. If these are now renamed and suitable changes are made to the 
length scales, a new point on the parametric may may be obtained. If a tilde 
is used to signify the conjugate run, we may write 

Y = 1/y (33) 


(32) 

at 

the 


*For rmlti.-ple-pa.-ix> oases^ see Section 5. 
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and 




(3h) 


The assumption implicit above, that true spiral form Is retained for the first 
half wave, has been confirmed in practice. However, attempts to carry out a 
similar procedure further down the wake have shown that true conjugates are 
not obtained; distortion of the spiral has become too significant. 


Size of time-step . The angle, A6, through which the joining line between 
two vortices rotates during one time step, At, is given approximately by: 


(Ti' + 12') 

A0 = — At radians 

2ir(nj - np)^ 


(35) 


For consistent accuracy, it is fairly obvious that the value of the stepping 
index, A0, must be maintained at an appropriately small level for all cases. 

It is also evident that as np or loading becomes large, a large number of 
small time steps will be required. Typical step sizes, At, In the previous 
study have ranged from 1/3 to 2.0, depending upon flap span, at a (Ci_/PR) 
value of 0.171^- It is generally found that the longer running cases (larger 
rip) also require the smaller values of time-step. This affects computation 
time severely.. 

Figure 4.13 shows the effect of step size on t ime-to-converge for four 
flap and tip vortex geometries. Quite good accuracy is obtained provided the 
stepping index is maintained at 8-degrees or below. Thereafter, the larger 
Hp cases degrade fairly rapidly. This occurs not only because longer times 
result at larger n^b'ut al-so because errors are larger at any given step size. 
Figure 4,l4, which shows the relative distance between two vortices after 50 
time units, demonstrates this. In the absence of center-plane images, an 
ideal stepping system would maintain constant spacing indefinitely. Figure 
4.14 shows that deviations from this trend increase with pp. 
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5. MULTIPLE-PAIR CALCULATIONS 


5.1 Introduction 

For each vortex pair in a multi -component system, there are at least 
three independent variables, even if a horizontal, unswept trailing edge is 
assumed. These are: vortex strength, span position, and core radius. Effec- 

tive viscosity is a possible fourth candidate. For N multiple vortices there 
are thus 3N or 4N minus one variables. (Length scale may be embedded in the 
core Reynolds number.) For a twin vortex system, a two variable map like 
Figure 4.7 must be repeated at least over a parametric range of Ci_/Al , for 
example, even for constant core properties. 

It is evident that it is not profitable to attempt to extend generalized 
results beyond twin-pair cases. A number of realistic, multiple-pair cases 
have therefore been selected, as examples for the present section, involving 
up to seven vortex pairs. These cases, defined in Figure 5*1 » will be 
described further in later sections. 

An additional motivation for treating cases individually concerns the 
calculation itself. For example, if a 4-pair rather than a two-pair case is 
considered, point density must be increased along each vortex, by a factor of 
two, giving four times as many points and a factor of I 6 on the number of 
velocity calculations required. Yet a further factor of two is-needed to 
properly relate the time step-size to vortex spacing*, yielding an overall 
computation Increase of 32-times, or 2^. Computer run time considerations are 
much more critical than core storage requirements, which vary only as the 
square of the number of vortex pairs. 

For these reasons, and because of the relative cost-per calculation, 
three-pair runs and up were carried out on Lockheed-Georgia computers only as 
necessary for program checkout, i.e. one or two time steps. The multiple pair 
runs reported in this section were made on the NASA-Ames CDC 7600 computer 
during a three-week period which included program checkout. Part of the 
development of the aircraft-fixed algorithms was also carried out at this 
time. This was sufficient to complete a number of test cases and to establish 
preferred methods and guidelines. Owing to the "one-shot” nature of the runs 
it was not possible, however, to produce completed calculations which may be 
regarded as totally definitive. The Ames results, plus subsequent studies,, 
indicate that the program does have the required capability, however. 


*Thta is mandatory for aircraft-fixed calculations ^ but some relief is 
possible in space-fixed cases. 
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5.2 General Approach 


The decay of a complex, multi-vortex wake Involves successive reductions 
in the number of vortices, either via convergence at the center plane (see 
Section 4) or by merging of two vortices Into one (see ref. 10). Figure 5-2 
shows how the length of the wake may be partitioned on this basis. In the 
present context, the terms “neai — "middle- or "fai — wake" relate only to 
the number of vortices present and have no specific counterparts in miles. 

Starting at the downstream end, we see that the "far" wake contains just 
one vortex pair. Since this has the potential for very long life under calm 
conditions, it is important to design the initial wake configuration in such 
a way that this remaining pair is a weak one. (We saw in Section 4 that, with 
suitable starting conditions, earlier convergence of stronger pairs is 
possible) . 

The "middle" v^ake involves a twin pair. This may be dealt with either 
using parametric plots or on a specific basis depending upon whether the 
geometry is near-spiral or has suffered noticeable distortion in the near 
wake*. In the latter case, space-fixed calculation methods may be used, as 
described In Section 4 . 

It Is the "near" wake which is of principal concern in this section. 
Considering briefly the simplest, three-vortex case, it Is noted that an em- 
bedded twin-pair wavelength will no longer characterize the motion. In most 
cases, this will prevent effective use of a space-fixed analysis, because the 
wavelength over which the motion is cyclic (assuming no convergences) is 
impract I ca 1 ly long. It follows that an aircraft-fixed analysis must be used. 

The dependence of computer run time upon the fifth power of the number 
of vortices (see Section 5.T) makes it highly desirable to eliminate from the 
calculation those vortices which converge early and to merge those which 
spiral together In close proximity. Figure 5.3 shows examples of both types. 
This particular case concerns a "sawtooth" loading configuration designed by 
Rossow (ref. 5 )* The reduction from seven to four vortex pairs is clearly 
most desirable, computationally. 

The results shown In Figure 5*3 were obtained using an aircraft-fixed 
two-dimensional calculation. A corresponding three-dimensional relaxation 
was attempted but proved unstable. Since complicated 3-D relaxation cases 
are likely to require individual treatment to overcome stability problems and 
are liable to be very expensive to run, it may be acceptable to rely on a 
simpler method for the early wake provided that the convergence or merging 
events are strongly indicated. If this is so, the predictive accuracy of the 
overall calculation will probably be acceptable. 


■'If speoif'i.Q acre details are to be used^ the parametric approach is likely 
to be precluded. 
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In summary, the overall approach includes the following steps, given 
vortex dispositions along the trailing edge and their strengths: 

(la) Use 2-D, reference 9, or (preferably, if stable) 3“D relaxation 
methods to identify and treat each vortex merging and convergence to the center 
plane. This will usually reduce the wake to no more than three or four vortex 
pairs. 

(lb) Continue to treat the remaining wake using aircraft-fixed methods 
until two pairs remain. 

(2) If the remaining twin pair is largely undistorted and individual 
core details are not available, a parametric map may be consulted, A recom- 
mendation for the preparation of such plots Is included later in this report. 

( 3 ) Once one member of the twin pair has either converged on the center 
plane or'merged with its mate, the ensuing single pair wake cannot be calcu- 
lated by present methods because it is straight and undisturbed in still air. 

If it is possible to postulate a periodic disturbance, the methods of Section 
3 may be used. 


5*3 Three-Dimensional Relaxation of the Near Wake 


•Before starting a three-dimensional relaxation, it is desirable to 
determine a reasonably close approximation to the final result. This reduces 
computer time and may prevent instability problems in some cases. Either a 
two-dimensional or a quas i -two-dimens ional (ref. 9) initial run may be used. 
Figure 5*^ shows an example of each type, together with a fully-relaxed, 
three dimensional run for the same case. Both of the two-dimensional methods 
display wake geometries very similar to the fully-relaxed, three-dimensional 
method. The fully-relaxed run does not extend as far downstream as the 
two-dimensional runs due to computer limitations. 

The three-dimensional relaxation process is well behaved for all the 
three-pair wakes reviewed. Disturbances originating near the trailing edge 
propagate away at approximately mainstream speed (see Figure 5*5) leaving a 
stable vortex system whose transverse velocities are of order 10"^ or 10”^ — 
the limit of computer resolution. The vertical arrows in Figure 5-5 depict 
the downstream end of the converged region. 

Attempts were made to converge 4-, 5~s 6-, and y-patr cases in the same 
manner as above, with decreasing degrees of success. However, all of these 
cases employed a swept trailing edge starting line which caused elements on 
adjacent vortices to be staggered relative to each other. It is believed 
that this destabilized the relaxation calculation. 
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5*4 Three~Pair Cases 


Certain three-pair wakes may be described as neai cyclic [see Figure 
5.6(a)], while in others [Figure 5.6(b)] the third vortex causes irregulari- 
ties in the shorter wavelength which gives the overall motion a more transient 
nature. in such cases, there is, therefore, no prospect for using a space- 
fixed analysis which uses reflection techniques based on the shortest- wave- 
1 ength. 

in neai — cyclic cases [Figure 5.6(a)], an attempt was made to simplify the 
problem by assuming that beyond the initial relaxation region, the third 
vortex remains straight/' A space-fixed calculation at the shorter wavelength 
was then carried out with the full streamwise image system, as used previously 
for twin pairs. As for the corresponding aircraft-fixed case, the innermost 
vortices travelled together as they were convected downward and outward by 
the main flap vortex. However, a rapid spiraling then occurred which is un- 
characteristic of the aircraft-fixed result. The straight-vortex approxima- 
tion was therefore unsuccessful. 

Recognizing the inherent aircraft-fixed nature of the expected phenomena, 
the above data were scanned as If by an observer moving at aircraft speed. 

This produced improvements in some cases, but cyclic instabilities which 
occurred at frame passage frequency prevented useful results from being ob- 
tained. However, deeper investigation suggested a hybrid aircraft-fixed/ 
space-fixed calculation technique which appears promising. This will be 
described in Section 5.6. 


5.5 "Jumbo-Jet‘' and "Sawtooth*' Cases 


The Boeing 7^7 runs outlined in Figure 5*1 were selected in the light 
of recent wind tunnel and flight wake hazard test experience at NASA concern- 
ing the effect of removing the outer flap and using the inner flap alone to 
restructure the vortex wake, and concerning the effect of landing gear deploy- 
ment on such results. In the present calculations the representation of the 
gear is via its estimated effect upon span load: this omits all viscous wake 
considerations. Nonetheless, the geai related span-load changes are found to 
have a significant impact upon the vortex wake structure. Pre-runs, to 
identify convergence, merging, etc. (see Section 5.2), permitted simplifica- 
tions to the distributions shov/n in the upper part of Figure 5-1 such that 
three-pair loading configurations could be used as depicted in the lower part 
of that figure. Three-dimensional relaxation of these three-pair configura- 
tions to a steady state then produced the wake structures shown in Figure 5-7. 

The upper two cases of Figure 5-7 correspond to early flight tests which 
were carried out, geat — up, for full flaps and for inboard flaps only. Wider 
vortex dispersion is evident In the latter case, which should be beneficial. 
Lowering the gear reduced the strength of the (anti-lift) Inboard flap vortex 
in both cases. With outer flaps present, lowering the gear caused a much 
less regular periodic near-wake solution. The reduction in inboard vortex 


^Though not fixed. 
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strength, due to the gear, also greatly reduced the tendency of the flap and 
tip vortices to spread apart. A more hazardous wake is therefore implied with 
the gear down. This corresponds to flight test observations. 

The results of Figure 5-7i and other '‘jumbo’'"jet runs not reported here 
because of step-size problems, have not been extended to determine time-to- 
converge. Two-dimensional extensions of the B7^7 runs have been made which 
indicate that convergence is occurring in the cases shown, but such runs can 
give no further information because curvature terms, which lead to convergence, 
are lacking. 

It is urged that the ‘'jumbo‘‘-jet demonstration runs started during the 
present study, are completed and compared with corresponding flight and wind 
tunnel test results. The use of the sliding-frame calculation method (see 
below) should make such studies more economical than was envisaged 
previously. 


5.6 The Sliding-Frame, Aircraft-Fixed Calculation Method 

In previous sections of this report, space-fixed and aircraft-fixed 
calculations have been regarded as quite se'parate approaches for solving the 
far vortex wake. However, a review of results, particularly the moving ob- 
server, space-fixed study of Section 5-^, has suggested a hybrid method which 
is best regarded from an aircraft-fixed point of view. 

In Figure 5-5, the history of an aircraft-fixed three-dimensional relaxa- 
tion is dominated by a disturbance wave which propagates downstream at approxi- 
mately mainstream speed and leaves behind the desired steady solution. Let us 
define the wave propagation speed as 3 times aircraft speed. We may choose a 
time, in the aircraft-fixed calculation, at which direct wing influence is no 
longer felt at the position of the disturbance front: "TIME =20" is suitable 

(Figure 5.5)- A calculation frame, stretching, say from Station I 6 to 
Station 28, is now drawn around the disturbance front. If this frame is slid 
along the wake at 3 times mainstream speed, the disturbance front in Figure 
5*5 remains centered and converged wake Is always contained In the left half 
of the frame. If suitable precautions are taken with conditions at the ends 
of the frame, particularly upstream, this provides an attractive and economical 
basis for a wake solution. In fact, if g = 1, the method reduces to a space- 
fixed calculation. The corresponding (3 = 1) program changes are trivial. 

Treatment of the downstream end of the frame presents little difficulty; 
a two-dimensional type of initial solution will clearly suffice, as for the 
aircraft fixed case. At the upstream end care must be taken both to ensure 
that a sufficient length of converged wake is included and to provide suffi- 
ciently realistic upstream conditions that the converged wake within the box 
is not contaminated. Such questions could be resolved by comparison with 
existing cases run with the current aircraft-fixed algorithms. 

Since the program changes involved are minor, and test examples exist, 
the above modifications and checkout appear to be highly worthwhile. 
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6. CONCLUSIONS AND RECOMMENDATIONS 


A computer program has been developed which uses finite-element three- 
dimensional relaxation methods to calculate the development of vortex wakes 
behind aircraft for a considerable downstream distance. The program is now 
operational at NASA-Ames. The inclusion of a self-induction term in the 
solution, dependent upon local curvature and vortex core radius, permits 
calculation of finite lifetimes for systems for which infinite life would be 
predicted two dimensionally. The present report describes the program to- 
gether with single-pair, twin-pair, and multiple-pair studies carried out 
using it. Though "jumbo'' aircraft wake examples are included, these have not 
been carried far enough for other than qualitative conclusions to be clrawn- 

The wake studies have lead to the follov/Ing conclusions; 

Si ngle-Pai r Wakes 

(1) There is a lower limit to the wavelengths at which the "Crow"-type 
instability can occur. Below this limit, self-induction effects cause the 
plane of the disturbance waves to rotate counter to the vortex direction. 

(2) Disturbance waves in a horizontal plane are several times more 
stable than the classic, A5“degree "Crow" case. Vertical and spiral waves 
are slightly more stable. 


Twi n-Pa i r Wakes 

(3) There is an important distinction between single-and twin-pair 
vortex convergence at the center plane. The first is usually a space-fixed 
event, but the latter moves at aircraft speed. 

(A) The self-induction term causes the time taken for vortices to con- 
verge at the center plane to be a strong function of wavelength (see Figure 
4.10). Either vortex can converge first (Figure 4.6). 

(5) in the limited twin-pair studies performed for constant C^/P^, 
t ime-to-converge increases with increasing flap span. 

(6) Choosing equal strengths for flap (for np ^ -333) and tip vortices 
maximizes the strength of the vortex which survives after one has converged 
at the center plane (see Figure 4.9)- 


Multiple-Pair Wakes 

(7) Though space-fixed axial-imaging techniques can be used for single- 
and twin-pair wakes, they are not applicable directly to multiple pairs. 
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Calculations must at least start in an aircraft-fixed coordinate system. 
Thereafter, a "si I ding-frame" method is suggested which has potential for 
applying three-dimensional relaxation methods economically to the far wake. 

The method is very similar to a space-fixed calculation but is controlled 
differently (see Section 5 . 6 ). 

(8) Pilot studies of Boeing 7^7 configurations show correct qualitative 
response to removal of the outer flap and to gear deployment, as compared with 
wind tunnel and flight test experience (see Figure 5 . 7 ). 


Recommenda t i ons 


( 1 ) More detailed studies should be undertaken of the separate effects 
of vortex core radius and relaxed wavelength on the time to attain center- 
plane convergence in twin-pair cases. This should provide useful design 
guidance for low-hazard wakes. 

(2) Further time-to-converge, twin-pair parametric studies (see Section 

k.2) should be carried out, holding constant for each of several maps 

like Figure 4.7a and 4.8. 

(3) The "sliding frame" calcultion method should be implemented and used 
to complete the present jumbo-jet studies. 

(4) More comprehensive studies should be made concerning wakes designed 

for early disruption, e.g. the "sawtooth" design. This requires more detailed 
aircraft fixed runs for near wake, including an improved rationale for vortex 
merging. If significant vortices then remain, the '^H^r^g J"£am^j]ie^ 

_shoujd*-be -appHed-;^ — - _____ 
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APPENDIX 


VORTEX LOOPING COMPUTER PROGRAM OPERATION 
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1. GENERAL PROGRAM OPERATION 


A computer program has been developed to apply the previous methods to 
multiple vortex pair instabilities. The program has been implemented on the 
UH06 computer and the MAC-16 mini-computer at Lockheed-Georg ia and the IBM 
360-67 and the CDC- 76 OO computers at NASA-Ames. 

The program has the capability to analyze up to seven vortex pairs and 
up to ninety points per wavelength in the cyclic, space-fixed option and ninety 
points per vortex in the aircraft-fixed option, in the space-fixed option 
each point is perturbed at each time step and a new geometry is calculated for 
the next time step. Since the frame of reference is fixed in space, no veloc- 
ity Is imparted to the system except that induced on each point by the system. 
The calculations are continued until the vortex system becomes unstable. 

The aircraft-fixed option is more complex. An initial two-dimensional 
time dependent, space-fixed calculation is performed. Time is then converted 
to downstream distance to set up the initial geometry for a three-dimensional 
time relaxation. In the aircraft-fixed frame of reference the free stream 
velocity must be added to the Induced velocities. The free stream velocity 
causes the points to move downstream. To keep the frame of reference fixed 
with the aircraft, the points must be reordered at each time step. Each 
point, N, at time, t, becomes the point, N+1 , at time, t+At. The points on 
the wing are fixed with time, and the trailing edge point is used twice at 
t+At,- once fixed on the wing and once shed downstream. The last downstream 
point is removed at each time step. The calculation continues until equi- 
librium is reached. This normally requires that the number of time steps 
equal the number of points per vortex considered. 


2. INPUT 


The general input parameters to the vortex looping program are; 

Spatial coordinates (x, y, z) 

Vortex strength (r) 

Semi -wave length (X/2) 

Wave amplitude (a) 

Core radius (6) 

A description of the Input parameters is given in Table A-1 . A1 l~var tables 
used in the vortex looping program are normalized on wing semispan (b/2) and 
free stream velocity (U^) • if the subscript, s, denotes normalized value 
used in the program and the subscript, o, denotes dimensional values, then; 
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TABLE A-1. INPUT PARAMETERS 


IPRT 

NCASE 

ITAPE 

lU 

IW 

IWR 

IPRTA 

IBACK 


Frequency of printout 

Number of cases to be run (l in current version) 
Graphics or plotter option (1 If desired) 

Output unit for graphics or plotter 
Output unit for printout 
Printout option (l if desired) 

Frequency of output for graphics or plotter 
Backward march option (1 if desired) 


J 

DELTMC J) 
RNCM< J) 
IPM( J) 
IZLMC J ) 

NPM( J ) 
MTIMM ( J ) 
ICASM( J ) 
lOPTMC J ) 


I0P2MC J ) 

TIMEM( J) 
DELTBC J ) 

CBAR( J) 
MPT5( J) 

TITLEdl. J> 


Case number 

Time increment 

Reference Reynolds number 

Number of wavelengths considered 

Number of points. per wavelength — must be even if 
I0P5M(J)=D. For aircraft-fixed or 2-D options 
IZLM(J)=1- 

Number of vortex pairs 

Number of time steps 

^aphics case number (not used in current version) 

Input option 

0 Sinusoidal input 

1 Point-by-point Input 
-1 Spiral input 

Full/half wave option 

0 Half wave considered 

1 Full wave considered 

Initial time 

Time increment (= distance increment) for initial 2-D 
aircraft fixed calculation. DELTBC J) must equal 
DELTMC J ) . 

Wing mean aerodynamic chord. 

Number of time steps for initial aircraft fixed calcula- 
tion. If MPTSCJ)=0, then space fixed option is used. 

Title 
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TABLE A-1. INPUT PARAMETERS (Concluded) 


I 

NPAMd » J) 
CIRCMCIr J) 
WVLMd . J) 
AMPXMdi J) 
AMPYMd , J) 
RCORMd f J ) 
SSMd. J), 
HINMdf J) 
PHASMd . J.) 

H 

XINCNrI, J) 
YIN<N,I, J) 
ZIN(N,I, J) 
RCOMKN.Ir J) 


Vortex pair number. For aircraft fixed calculations, 
vortices must be input from inboard to outboard. 

Vortex reference number 

Vortex strength 

Wavelength/2.0. Must be the same for all vortices. 
Amplitude of slnusoldal/spi ral disturbance in x-dlrection 
Amplitude of sinusoidal/sprial disturbance in y-direction 
Core radius •- 

Initial spanwise location of vortex 

Initial vertical location of vortex 

Phase shift of disturbance. Must be the same for all 
vortices. 

Point number for point-by-point input 
X location of N 
Y location of N 
Z location of N 
Core radius of N 
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>■5 = V(U.b/ 2 ) 

V: = V; /(U ) 

' S ' O “ 

ts = to(uy^/2) 

>^s = V(b/2) 

as - ao/(b/2) 

= S^/{h/2) 

The input format is summarized in Table A-2. Note that x is the spanwise 
coordinate, y is the vertical coordinate, and z is the downstream coordinate. 

3. MAIN PROGRAM/SUBROUTINE FUNCTIONS 

The principal functions performed by the main program and each sub- 
rout i ne are: 


SUBROUTINE 


FUNCTIONS 


Main Program 

1 . 

Case data input. 



2. 

Initial geometry for one side of each 
determined using sinusoidal or spiral 
or point-by-point input. 

vortex pai r 
distribution 


3 . 

Other parameters used in calculations 
(XjDUM, CXj=rAxj, etc.). 

determined 


k. 

Displacements and new geometry for each time step 
calculated. 


5. 

Geometry at each time step output. 


POOR 

1 . 

Induced velocities summed. 



2. 

Element compression applied. 
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TABLE A-2, INPUT FORMAT 


Card Format Input Parameters 

1 Rl'f IPRT, NCASEr ITAPE, lU, IW, IWR. IPRTA, 

IBACK 


2 SFlO.OfTIH- DELTM(J), RNCM(J), IPMCJ), IZLM(J), NPM<J) 

MTIMM(J)r ICASM(J), lOPTMCJ), lOPSM(J) 


3 3F10.Q.I4- TIMEM(J)r DELTB(J), CBAR<J)r MPT5(J) 

(If aircraft fixed oaZaulations are not to be done^ 
then MPTS(J) =0.) 


h EOAif TITLEdIfJ) 


If I0PTM(J)=0 


5 It^.SFR.D NPAM(IrJ), CIRCMCIiJ), WVLM(I,J). AMPXM<I,J), 

AMPYM('l,J)r RCORMCI.J). SSM<I,J). PHASMd.J). 
HINMd . J) 

~ — — (Repeat~card 5~fdr eddhTvdrtex pdir7)~~ ~ ~ 


If IOPTM(J)=l 


5. I^frSFS.O NPAMd.J), CIRCMdrJ). WVLMCIrJ) 

6 ‘fFlO.O XINCN.IrJ), YIN<N,IrJ), ZIN(N,I,J), 

RCOMKN ,I , J) 

(Repeat card 8 for each point (N~lj INPI) on a 
single vortex pair (one side only).) 

(Note: Calculations are performed on the points 

= INPI-1.) 


(Repeat cards 5 and 6 for each vortex pair. ) 


Sh (Repeat cards 2 through 6 for each case.) 



FAAVEL 


CURV 


STAT 


A list 
f o 1 1 ows . 


1, Induced velocities from each point determined. 

2. imaging techniques used to determine the effects of up 
and downstream points and the vortices shed from the 
other v/ing. 

1. Local radius of curvature calculated. 

2. Core radius at time, t, calculated. 

3. Self-induced velocity calculated. 

1. Final geometry at end of backward march compared to 
initial geometry. 

4. PROGRAM LISTING 

i ng of the vortex looping program as used on the U1106 computer 
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1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 
-34 - 

35 

36 

37 

38 

39 

40 

41 

42 
,43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 


C PROGR.AM MAIN ( IN^UTf OUTPUT? TA'^E5=INPUTf TA»E6=0UTPUTr T.A'^ES t T.APE1 ) 

C MAIN PROG 

DIMENSION NPAM{7?1 )?CIRCM(7?i )?WVLM(7»1 )?AMPXM(7?1 )? 

1 ,AMPYM(7?1 )?RC0RM(7?I )fSSM(7?l )?PHASM(7?1 )?ICASM(1 )r 

2 MTIMMd )fDELTM(i )?RNCM(1 )?I°M(1 )rI2LM{l ),NPM{1 ) 

3 ?XIN{92r7»l )?YIN(92?7?1 ) ? ZIN,( 92? 7? 1 5 ? RCOMI ( 92? 7 ? 1 ) 

4 ? ioPTMCl) , fopEMd) ?HINM(7?i) ? TIMEM ( 1 ) ? TITLE(20 ? 1 ) ?TTL{ 20 ) 

5 ?D£LT8(1) i|f^PtS(l) ?CqAR(l) . 

^DIMENSION X«(92?7) ?YP(92?7) ?ZP(92f7) ? NP.A ( 7) ? WVL ( 7) ?,AMPX(7) ?.AM°Y(7) 
1 ?RC0RE(7) ?SS(7) ?SEMIW(7) ?=>H(7) ?°HAS(7) ?HIN(7) 

.COMMON, X( 92? 7) ^ Y( 92? 7) ? Z{ 92?7) ? U( 92): 7) ? V ( 92? 7) ? W ( 92?7) ? , 

1 XDUM{93?7) ?YDUM(93?7) ?ZDUM(93?7) ? CX ( 92? 7) ? CY-(92 ?,7) ? CZ ( 92 ? 7) ? 

2 CIRC(7) ?,ALMBDA{7) ?x6(92?7> ?Y0'{92?7) ?Z0(92?7) ?DELRM{92?7) 

3 ?,RC0M(92?7) . , . , - . . , . 

DATA “I/3.1415926/? QPI/.7853982/? Q°I2/1. 5707963/? ZERO/0. OX ?PI2/6*2 
$831852/ 

90 CONTINUE 
IR=5 
I2D=:6 

C IN°UT „ - . , . . 

READC IR? 106) IPRT?NCASE? ITA°E? lU? IW? IWR? IPRTA? IBACK 
IF(IPRT.LT.O) |30 TO 998 
DO 54.J=1?NCA5£ 

R£AD(iR?107)D£LTM{ J) ?RNCM{J) ?IPM( J) ? fZLM{ J) ? NPM( J) ?MTIM'1(J) ? 

1 ICASr^(J) ?IO»XM(J) ?I0P2M{J) 

READ (IR? 120) TIMEM (J) ?DELTB( J) ?CBAR( J) ?M«TS( J) 

R£AD( IR?302) (TITLEdl? J) ? 11=1? 20) 

NPRE=NPM(JJ 

IF( IZLM{ J) .LE.i) I2D=1 
DO 54 I=1?NPRE 
IOPTT=IOPTM( J) 

IF(IOPTT.LE.O) GO TO 56 

READ ( I R ? 1 0 3 ) N PAMI 1? U )-?-C I RCM ( I ? J-) ? WVtM(IrJ) 

INPi=IZLM{ J)/2+2 ' ■ 

iFCI0P2M(J) .GE.1) INPI=IZLM{ J)+2 
iF(I2D«G£.l) INPI=1 
DO 57 K=l?INPi 

57 READ { IR ? 117 ) XIN ( K ? I ? J ) ? YIN C K ? I ? J ) ? ZIN ( K? I ? J ) , RCOMI { K ? I , J) 

GO TO 54 

56 .READ(IR?103)NP,AM(i ? J) ? CIRCM ( I ? J) ? WVUM ( i ? J) ? AMPXM { I ? J) t 
. 1. AM°YM(i?J) ?RCORM{I?J) ?SSM(I?U) ?‘^HASM(I?J) ?HINM(I?J) 

54 CONTINUE 

DO 5 J=1?NCAS£ 

C INITIAL .COND, FOR EACH C.ASE 

DO 83 1=1? 20 ' ■ 

88 TTL(I)=TITLE( I? J) 

D£LTZ=0.0 
I2D=0 

MTIM£=MTIMM{ J) 

IAC=0 
JAC=0 

IF(MPTS(J) .GE.l) iAC=i 
IFdAC.GE.l) JAC=1 
IFdAC.GE.l) MTIME=MPTS(U) 

XMAC=CBAR( J) 

ICASE=ICASM( J) 

D£LT=D£LTM{ J) 

IF{ lAC.GE.l) DELT=DELTB( J) 
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59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 
75 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 
87 
83 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 
117 


RNC=RNCM( J) 

IP=I'=MC J) 

6 IZUINK=IZLM( J) 

IF(IZLINK.tEci) I2D=1 
NOAIR=NPM( J) 

I0PT=I0PTM( J) 

I0°2=I0P2M{ J) 

DO 55 I=1»N=AIR 
NPA(I)=:N°AM( If J) 

IFdOPT.UE.O) GO TO 58 
INP=IZLINK/2+2 
IF(I0P2.6E.l) IN‘==iZLINK+2 
IF(I2D.GE,1) 

DO 97 K=1»INP 
X(Kd)=XIN(Kf IfvJ) 

Y(Kf I)=YIN(Kdr J) 

ZCKf I)=ZIN(Kf . 

ROOM ( K f I) =RCOMI ( K f I f J ) 

X0{Kr I)=XIN(Kfif J) 

YOtKf I)t:YIN{kdf J) 

97 ZOCKf I)=ZIN(Kf If J) 

58 CONTINUE 
WVL(I)=WVLM(if J) 

IFCIOPT.GT.O) GO TO 55 
AMPX( i)=AMPXM( If J) 

AMPYC I)=AMPYM( If J) 

RCORE( I)=RCORM('if J) 

S3(I)=SSM(If J) 

PHAS( I}=PHASM(Jf 
HIN( i)=HINM( If J) 

55 CIRC{i)=CIRCM( I/J) 

RNC=RNC+10000C:. 

IMARCH=0 

iF{iWR.L£.0)GO TO 61 
WRITECiWf 303) 

WRITEdWf 302)TTL 
WRITECIWflOS) 

IF(D£LT}999?999f66 . , 

60 V/RITEdWf 109)DELTfRNCf i='f IZLINKfNPAIR 
IFdAC.LE.O) 60 TO 26 
WRITECiWf305)XMAC , 

26 IFdOPT.LE.O) GOdO 67 
WRITEdWf 118) 

GO TO 68 

67 WRITEdWf 114) 

68 CONTINUE 

DO 52 I=lfN°AIR 
IFdOPT.LE.O) 60 TO 59 
WRITE( IWfl04-)IfCIRC(i) 

GO TO 52 

59 WRITEdWf 104) I f CIRC (I) f WVL( I) f AN^’Xd) f AMPYCI) f RCOREC 
1 fSSd)f°HASd) 

52 CONTINUE 

61 RIZL=IZUNK 
IT£=0 

IST=0 . „ „ 

C INITIAL, VORTEX GEOMETRY ORIGMAE PAGE K 

DO 66 I=lfNPAiR OF POOR QUAIOT’ 

IF dOP2)70f70f71 


37 


‘HI 



118 

119 

120 
121 
122 
123 
124- 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 
161 
162 

163 

164 

165 

166 

167 

168 

169 

170 

171 

172 

173 

174 

175 

176 


71 IHALF=IZLINK 
GO TO 72 .. 

70 IHALF=izU INK/2 

72 IFUHALF-gO) iO^lO^ll 
11 .WRITE( IWfllO) IZLINK 

GO TO 999 

io IHALFP=,IHALF+2 . 

IF(I2D.GE.l) IHALFPrl 
IHPl^zIhlAEFP+l 
IH1=IHALF+1 
rBR:=(2+IHl)/2 
IF(I2D.GE.l) GO to 82 
IFCIOPT.LE.O) GO TO 63 
DO, 64 N=2riHALFP 
NN=N-1 

XDUN(Nf I]=0.50+(X(N» I)+X(NM» I) ) 

YDUN(Nf I)=d,50*(Y(Nr I)+,Y{NM» I) ) 

64 ZDUMCNf i)=0.50+(Z(Nf I)+Z(NM»I) ) 

DO 65 N=2rIHl 

NP=N+1 

■ CX(Nf I)=CIRC( i)*(XDUM(N» I)~XDUM(N°r I) ) 

CY(N»I)=CIRC( I)*(YDUM(Nf I)-YDUM(NP^ I) ) 

65 CZ{Nf i)=CIRC{Xi*(ZDUM(N» I)-ZDUN(NP/l) ) 

CX{lrI)=-CXC2r I) 

CY{lri)=-CY(2tI) 

CZ(lrI)=CZ{2»I) 

CXC IHAUFP» I)=-CX( IH1» I) 

CYC IHALFPf I)=-CY{ IHlf I) 

CZ ( IHALFP r i; =CZ ( IHl r I ) 

ALMBD A ( I ) =2 . 0 * { ZDUN ( IHALFP f I ) -ZDUM ( 2 » I ) ) 

IF { I0P2 . GE . 1 ) ALM8DA ( I ) =ALNBD A ( I ) /2 . 0 
GO TO 66 
63 AN°=AM‘^Y(I) 

ALMBDA { I ) =2 . 0’f=WVL ( I) 

_ .P_Ha-)-=pHAS(jl*AL-MBDA{-ip ~ 

ZLI NK= ALMBD A ( I ) /R I ZL 
HZLINK=ZLINK/2.0 
SEMIWC I)=SS(i) 

ZDUM C 1 r 1 3 =~ ALMBD A ( I ) /2. O+PH { 1 3-ZLINK 
IF(IO'=T.,LT.O) GO TO IS 

YDUM ( 1 r I ) =HIN ( I ) +AMP+COS ( PI2’»=ZDUM ( 1 r I ) / ALM8DA ( I ) ) 

GO TO 19 

18 YDUM ( i r I ) =HIN ( I ) +AMP+SIN ( PI2*ZDUM C 1 f I ) / ALM8DA { I ) ) 

19 XDUM( 1 » I ) =SEMIW ( I ) +AMPX ( I ) +C05('^I2't=ZDUM ( 1 r I ) /ALMBDA ( I ) ) 
DO 12 N=2rIHPP 

NM=N-1 

ZDUMCNp I)=ZDUM(NM f i)+ZLINK 
IF(10°T.LT.Q) GO TO 27 

YDUM ( N > I ) =HIN { I ) +AM'3+C0S ( PI2*ZDUM ( N f I ) /ALMBDA ( I ) ) 

GO TO 28 , 

27 YDUM ( N? I ) =HIN { i ) +AMP+SIN ( °i2*ZDUM ( Nr I ) / ALMBDA ( .1 ) ) 

.28 XDUM(N» I)=S£MIW{ I)+AMPX{I)*COS(PI24=ZDUM(Np I)/ALM8DA(I) ) 
Z(NMf I)=0.50*{ZDUM(NMf I)+ZDUM{Nf I) ) 

IF(iOPT.LT.O) GO TO 29 . , 

YCNMf I)=HIN(I) +AMPtC0SC'-^I2'*=Z(NM>I)/ALM8DA( I) ) 

'GO TO 36 

29 YCNMf I)=HIN(I) +AMPmsSIN(“I2’«=Z(NM^I)/ALMBDA( I) ) 

36 XCNMf I)=SEMIW( I)+AM«XCI)*C0S{'^I2*Z{MM»I)/ALMBDA( I) ) 
XOCNMf D^XCNMr I) 
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177 

178 

179 

180 
181 
182 

183 

184 

185 

186 
187 
186 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 

199 

200 
201 
202 

203 

204 

205 

206 

207 

208 

209 

210 
211 
212 

213 

214 

215 

216 

217 

218 

219 

220 
221 
222 

223 

224 

225 
225 

227 

228 

229 

230 

231 

232 

233 

234 

235 


Y0{NMri)=Y(NM,i) 

ZOCNMf I)=2(NMf 1} 

CX(NM» I)=CIRC(i) ’H(XDUM{NMr ) 

CY ( NMr I ) =CIRC ( I ) ■*= ( YDUM ( NMf I) -YDUM ( Nr t) ) 

12 CZ(NMf I)=CIRC{I)>f=(2DUM{NMri)“2DUM(NrI) ) 

GO TO 65 
82 CX(lrI)=0. 

CYdr I)=0. 

CZCir I)=-CIRC(I ) 

ALMBDA(I)=2.0=*=WVL(I) 

66 CONTINUE 

IF(I2D.G£,1) GO TO 87 

IF(10'=>2.GE.l) D£LT2=(Z(IH1»1)+Z(2rl) )/2.0 
D£LTZ0=DELTZ 
87 CONTINUE 

WRITER IWrlll) 

I°TI=0 

IRTA=0. 

TIME=TIMEN{ J) 

DO 25 I=lrNPAIR 
' DO 25 N=lrIHALFP 
XP{NrI)=X{Nri) 

YP(NrI)=Y(NrI) 

ZP>(N> I)=Z(NrI) 
iF(IOPT.GT.O) GO TO 25 
RCOM(Nr I)=:RCOR£(i) 

25 D£LRM(Nf I)=0.0 
XBAR=0 . 0 
YBAR=0.0. 

deltyb=o . 6 

23 CONTINUE . 

IF(IAC.LE.O), ,G0 TO 2 
FAC=0.75*XNAC/'2.0 
FAC1=FAC^ ' . 

DO 80 I=lfN'=AIR 
KK=1 

RCOMldrlf J)=RCOM{lrI) 

YIN(1» If J)=Yb(lf I) 

ZINdf If J)=Zb(lf I)-2.0=J=FAC1 
IFd.EQd) , GO T0..81 
XlNdf If J)=X0df I-l) 

GO TO 80 

81 XINdf If J)=b.O - 

80 WRITE(IWfll2)KKfXIN(lf if J) f YIN { 1 f I f J) f ZIN( 1 f I f J) 
DO 3 I=:ifNPAIR 
DO 3^K=2f3 

RCOMI (Kf if J}=RCOMdf i) 

XIN(Kf if J)=X0(If I) 

YlNCKf If J)=Y0df I) 

ZIN(Kf If J)=Z0df I)-FAC~FAC1 

WRiT£(IWfll2)Kf XIN(Kf i f J) f YINCKf If J) f ZIN(Kf If J) 

3 FAC=“FAC 
2 CONTINUE 

DO 50 IT=lfMTINE 
C OUTPUT OF WAVE GEOMETRY 

IFdAC.GE.l) K=IT+3 
IPCO— 0 

IFdMARCH.GE.DGO TO 30 
IF( IWR.LE.O)GO TO 20 


OBIGBIAS mgww . 

POOR QUAU:^ 
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236 

237 

238 

239 

240 

241 

242 

243 

244 

245 

246 

247 

248 

249 

250 

251 

252 

253 

254 

255 

256 

257 

258 

259 

260 
261 
262 

263 

264 

265 

266 

267 

268 
269 

-2-70 - - 

271 

272 

273 

274 

275 

276 

277 

278 

279 

280 
281 
282 

283 

284 

285 

286 

287 

288 

289 

290 

291 

292 

293 

294 


IF(IPTI.GE.IPRT) I‘^TI=0 

XBAR0=X8AR 

Y8ARb=Y8AR 

GXT=0,0 

GYt^D.d 

Gf=0.d^, 

DO 79 I=lrNPAIR 
DO 79 NM=2rIHi 
N=NN .... 

iF{I2D.G£.i) N=N?,!-i- ,, 

GXTsGXT+CZjNrl) ^=X(N» I) 

GYT=G YT+C2.( » I ) *Y ( N » I ) 

79 GT=GT+CZ{N»I) 

X8AR=GXX/GT 

YBAR=:GYT/GT 

VDR={YBAR-YBARQ)/DELT , 

PELTYB=DELtYB+VDR*DEL.t 
IF(IAC.GE.l) GO TO 7 
WRITECIWfll3)TlME 

7 IF(IPTI.GE*l),,GO .TO 20 
‘ rFdAC.GE.l) GO TO 8 

WRI T£ C I W r 304 ) XB AR t Y8AR r VDR > DELTY8 

8 CONTINUE 

DO 78 i=l^N»AIR 
DO, 78 N^lrlHALFP 
U{Nri) = (X(Nf Il-XR(Nd) )/D£LT 
V{Nf I) = (Y(Nd)-YP(Nf I) )/D£LT 
78 W(Nf I)=(Z(N»I)-ZP(NrI) )/D£LT 
DO 51 Izl/NPAIR 
iF(i2D.G£.l) ZDUM(lfI)=b. 
iFilAC.GE.l) GO TO 9 
WRlt£(IWri0S) If DELTZ 

9 CONTINUE 

DO 51 NM=lrIHAUF 

N=NN+1 - 

IF{I2D.G£.l) N=NN 
IFdAC.LE.O) .GO TO 13 
IF{IT.UE.i) GO TO 51 
KP=K-1 

.WRITEdWf 112),KPf XINCKPf If J) fYINCKPf If J) fZIN(KPf If J) , 

1 U(NfI)fV{Nf I)fW(Nf I) 

GOdO 51 

13 WRITECiWf 112)NNfX(N, i) f Y(Nf I ) f Z ( Nf I ) f U (Nf I ) f V ( Nf I ) f W ( Nf I ) 
51 CONTINUE 
20 IPTI=1°TI+1 

IFdTAPE.L£.0),GO TO 30 
IF(I2D.GE.l) GO TO 30 
IFdT£.dE.l)GO TO 17 
DO 41 I=lfNf=»AIR 
DO 41 NN=lfIHALF 
N=NN+1 
IXN=X(Nf i) 
lYN=Y(Nf I) 
l2N=Z(Nf i) 

IF( lABSdXN) .GE.lOO ) GO TO 42 
IFCIABSdYN) .G£,i00 ) GO TO 42 
IFCIABSdZN) .GE.lOO ) GO TO 42 
GO TO 41 
42 IT£=1 
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295 

296 

297 

298 

299 

300 

301 

302 

303 

304 

305 
,306 

307 

306 

309 

310 

311 

312 

313 

314 

315 

316 

317 

318 

319 

320 

321 

322 

323 
' '324 

325 

326 

327 

328 

329 

330 

331 

332 

333 

334 

335 
356 

337 

338 

339 

340 
. 341 

342 

343 

344 

345 

346 

347 

348 

349 

350 

351 

352 

353 


C 


WRITE { I Wf 115) 

GO TO„17 
41 CONTINUE 

IRC IPTA.GE. I«RTA) I=TA=0 
IF(roTA.6E.i)G0 TO 17 
IFdT.GT.l) GO TO 86 
IF(TIME.GT.O.Ol) GO TO 17 
86 CONTINUE 

DO 15 i=lfN“AIR 
DO 15 NNrlflHALF 
N=NN-H 

YTAP=Y(Nf D-DELTYB 
IF(JAC.GE.i) YTAP=Y(Nfi) 

ZTAP=2{N» I)-D£LTZO 
15 WRITE(IUr300)X(Nf I) f YTA“rZTA« 

17 IPTArli^TA+l 
30 IF(IT.GE.MTIME) GO TO 50 
DO 43 I=1^NP,^IR 
DO 43 NN=lrIHALF 
N=NN+1 

' IF(I2D«GE.l) N=NN 

IF( ABStXCNf I) )-l00G0Q0. )44-»45r45 

44 IF( ABS(Y(N> I) )-10300D0. )46» 45^45 
46 IF(A8SCZ(Nf I) >-1000000. >43r45f45 

45 IST=1 
WRITE(IWfll6) 

GO TO 48 

43 CONTINUE _ 

'CALCULATION ,OF INDUCED VELOCITIES 
DO 40 l=lrN°AIR ’ - - ~ 

DO 40 NN=2riHl 

NC=NN 

N=NN 

' IF(I2D«GE,i) , , . . , 

CALL PCOR(Nf Ir IPCOf iPf IBRrNPAIRf IHl » IHALFP r DELTt IWf istf I0P2r 
i;jACfX«AC) 

IFCiST.GE.l) GO TO 48 
IF(I2D.GE.l) GO TO 33 
IFCJAC.GE.l) W(Nf I)=W(NrI)+l,0 
IF(JAC.LE.O) GO TO 69 
IF(N,LE.3) GO TO 83 

69 CONTINUE . , . . 

CALL CURVCNCfDELTf IHALFflPCO^ ITf RNCf If UCr VCf WC? IWr ist) 
IF(IST.GE.i) GO TO 48 ' ‘ - - • 

GO, TO 84 

83 UC=0. 

VC-0. 

wc=o. 

84 U(Nf i)=U(Nf D+UC 
V{Nf I)=V(Nf D+VC 
WCNf I)=W(Nf D+WC 

40 CONTINUE 

MOD. OF "INDUCED VEL. SO TH.AT VEL ALONG VORTEX = 0 
iF( JAC.LE.d) 'GO TO 96 
IFdAC.GE.l) GO TO 96 
DO 1 I=ifN»AIR 


DO 1 N=lfIHl 
IF(N.LE.2) GO TO 93 
RX=XDUM(N+lr I)-XDUM(Nf I) 
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354 


R Y=YDUM { N+1 f I ) -YDUM ( N r I ) 

355 


RZ=ZDUM ( N+i , I ) -ZDUM ( Nr I ) 

356 


RPSQ=RY*RY+RZ+RZ 

357 


Ro=SQRT(R°SQ) 

358 


RSQ=R°SQ+RX+RX 

359 


R=SQRT(RSQ) 

360 


AL31=RX/R 

361 


AL32=RY/R 

362 


AL33=RZ/R 

363 


AL21=0.0 ■ 

364 


AL22=RZ/RP 

365 


AL23=-RY/RP 

366 


SIGN=RZXABS(RZ) 

367 


ALil=SIGN*RP/R 

368 


ALl2=-siGN+(RX+RY)/,(R+RR) 

369 


ALi3=“SIGN* (RX*RZ),/(R+=RP)., . ^ 

370 


UP=ALll't=U ( Nr i ) +AL12<=V f Nr I ) -fALlS^f^W ( Nr I) 

371 


V'3=AL21+U { N r r ) +.Vl 22+ V ( f^l r I ) +/\L23* W ( N r I ) 

372 


WP=AL31*U ( N r I ) +AL32* V { N r I ) +AL33+^ W ( N r I ) 

373 


WP»’=W« . . . 

374 


U C Nr i ) =ALli +Uf=+AL21*V°+AL31+Wpp 

375 


V ( N r J ) =AL12 +UP+ AL22* VP+ AL32+ Wp° 

376 


W C Nr I ) =AL13*UP+AL23^VP+AL33+W^P 

377 


IF(N.GE.3) GO to i 

378 

93 

UCNr J)=0.b 

379 


V(Nr n=6.o 

380 


W(Nr,I)=0.0 

381 

1 

CONTINUE 

382 

96 

CONTINUE 

383 


DO 53 I=lfNPA,IR 

384 


IF{I2D.G£.l) GO TO 75 

385 


iF{I0P2)73r73r74 

386 

73 

Udr i)=U(2r I) 

387 


V(lrI)=V{2r I) 

_.388 

. 

- W ( 1-r IO-=”W (-2 r-iO“ - ~ 

389 


U(IHALFPri)=U(IHlrI) 

390 


V(IHALFP>rI)=V(IHlrIl 

391 


W(IHALF?rI)=-W{IHlrI) 

392 


GO TO 75 

393 

74 

iF(JAC,G£.l) GO TO 95 

394 


U(lri)=U(IHlri) 

395 


V(lrI)=V(IHlrI) 

396 


W(lr I)=W{IHlr I) ^ 

397 


U(XHAUFPri)=U{2rI) 

398 


V(IHALFPri)=V(2rr) 

399 


W(fHALFPrr)=W(2rI) 

400 


GO^to 75 . . , . . ... 

401 

95 

UCIHALFPr I)=2.0*U(IHlrI)-U(IHl-lr I) 

402 


V( IHALF'^r i) =2.0*V( IHlr I)-V{ IH^-lr I) 

403 


W ( IHALFP r I ) =2. 0 + W ( IHl r I ) -W { IHl-l r I ) 

404 

75 

CONTINUE . . _ 

405 


IFd.EQ.l) IPCO=IPCO-H , 

406 • - C 


CALCULATION OF VORTEX GEOMETRY 

407 . ' 


iF(IPC0-l)32r32r33 

408 • 

32 

DO 3i:N=lrIHALFP 

409 


XP'(Nr I)=X(NrI) 

410 


YP(Nri)=Y(NrI) 

411 


ZP(Nr I)=Z(Nr I) 

412 


X(Nr I)=Xa(Nr I)+U(Nr D+DELT 

42 





413 

414 

415 

416 

417 

418 

419 

420 

421 

422 

423 

424 

425 

426 

427 

428 

429 

430 

431 

432 

433 

434 

435 

436 
457- 

438 

439 

440 

441 

442 

443 

444 

445 

446 

447 

448 

449 

450 

451 

452 

453 

454 

455 

456 

457 

458 

459 

460 

461 

462 

463 

464 

465 

466 

467 

468 

469 

470 

471 


Y ( N f I) =YP ( N f i ) +V ( N f J ) ■4^DELT 
31 ZCNfi)=ZP(N^r)-t-W(NrD’4=DELT 

GO TO 

33 DO 34 N=iriHAUFP 

X(Nf I)=0.5f=(XPCNr|)-i-X{N> I)+U(N» i) ’t=DELT) 

Y(Nr I)=0.5’t=(Y°(Nri)-i'Y(Nf I)+V(Nr D+DELT) 

34 Z(Nf I)=0.5’t=(ZP(N> D-t-ZfNf I) +W(Nr I)*DEUT) 
iF(JAC.L£.6) GO TO 4 
IFdAC.GE.l) .GO 70 38 
DO 14 'N=lr IHALFP 
NDE=IHALF?-N-l-l 
IF( IMARCH.GE.l) NDE=N 

if(nde.i:e. 2 ) go to 4 

IDP=1 

IFdMARCH.GE.l) IDP=-1 
X C NDE-dDP d ) =X { NDE f I ) 

Y(NDE+IDPd)=Y(NDEd) 

14 Z(NDE-dD»f I)=Z(NDEd) 

4 CONtiNUE 

IFCJAC.LE.6) GO TO 38 
- DO 39 'N=l»3 
X{NrI)=X°(Nd) 

Y(N» |)=YP(Nd) 

39 Z(NrI)=Z°(Nd) 

IFdMARCH.LE.O) .GO JO 38 
X (I HALFP d ) =XS^ dHALRP d ) 

Y ( IHALFP r i ) =YP ( IHALFP r I ) 

ZdHALFPf I)=Z'^dHALFPd) 

38 CONTINUE 

IF(IAC.LE,b), GO to 35 
IFdMARCH.GE.l) .GO TO 35 
XIN(Kfi> J)=X{lrI) 

YIN(KrIfJ)=Y(lfI) 

ZIN(Kdf J}=ZIN{K-lf If J)4D£LT 
RCOMKKf if J)=RCOM(lf i) 

35 CONTINUE . ' ' 

IFCI2D.GE.1) GO TO 53 
DELTZ=(ZdHlf I)-i-Z(2fIJ )/2.0 
IFdOP2.LT.l) D£LTZ=b.Q 
IFCJAC.UE.O) Z(lr I)=Z{ld)“D£LTZ 
DO 47 'N=2f IHALFP 
NM=N-'l 

IF( JAC.LE.O) Z(Nf i)=Z(Nf D-DELTZ 
XDUM{NfI)=0.50*(X(Nri)+X(NMf I) ) 

YDUNUNf I ) =b . 50* { y'{ Nf i ) +Y( NMd) ) 

47 ZDUM{Nf I)=0.50*(Z(NrI)+Z(NMf I) ) 

DO 49 N=2dHl 

NP=N+1 . . 

' CX{Nd)=ciRC{i)*(XDUM(Nf I ) -XDUM ( NPr I ) ) 
CY(NfI)=CIRCd)*{ YDUM ( N r I ) - YDUM ( NP f I ) ) 

49 CZ(Nd)=CIRCd)*(ZDUM(Nf i)-ZDUM(N°f I) ) 

53 CONTINUE' . ’ , 

IFd^CO.LE.DGO to 30 
TIM£:=TiM£-KDELT 

50 CONTINUE , 

IFdAC.GE.l) GO TO 92 
DO 85 i=lfNPAIR 
DO 35 NpldHALFP , . , 

85 WRITE dWf 119) I»NfXCNfI)fY{Nd)fZ(NfI) f RCOMCNf I) 



472 

92 

CONTINUE 

473 

48 

CONTINUE , 

474 

99 

iF(,IWR.L£.0) .GO TO 98 

475 


WRITE(iWr30i) 

476 

98 

|F{IBACK.L£.6)G0 to 21 

477 


IF(iST.GE,3,) GO, TO 21 

478 


iF(IMARCH.LE.O)GO TO 22 

479 


WRITECIWf 113) TIME 

480 


WRIT£(IWf202) . ■ 

481 


DO 24. I=ifNPAIR 

482 


^am= 6 . b 

483 


XGMr'i . 0 

484 


WRlt£(IWf2b3 ) 1 

485 


DO 24 N=2fIHi. 

486 

24 

CALL StAT(Nf if IHlfXAMfXGl 

487 


GO TO 21 

488 

22 

IMARCH=1 

489 


D£LT=-D£LT 

490 


iPTI=IPTI-l 

491 


GO TO 23 

492 

21' 

IF(IAC.L£.6) GO to 5 

493 


IAC=0 ■ 

494 


iZLM( J)=MTIM£ 

495 


I2D=6 

496 


MTIM£=MTIMM( J) 

497 


DELT=DELTM( J) 

498 


GO TO 6 

499 

5 

CONTINUE 

500 

999 

GO TO 90 

501 

loi 

FORMAT ( 2b A4) 


102 FORMAT (2 13),, , 

103 FORMAt(i4»8F9.0) 

104 F0RMAT(I3»7Eil.4) , , 

105 FORMAT{iX»7HIPAIR =^l4»J^lm£LTZ„=f_Eil-A) 

1-0 6 'FORMAT C9 1 4')T~^ 

107 FORMAt(2F10.0r7I4) , . 

108 FORMAT(/flOXr4HD£LJ^ilXf5HRNREFf9Xr2HIPf3X/3Hi2Lf2?<»5HN'^AiR) 

109 F0RMAT(3X»2{2X»ei4.6),»3{2Xr,i4) ) 

110 FORMATC5(/X^829**+DiAQNOSTiC++i‘-T; — IZLINK MUST BE LESS THAN OR EQU 
SAL TO. 40. THE INPUT VALUE WAS (>I4>1H)) 

111 FORMAT(3(/) f5X»lHN»4X»4HX{N) »8X>4HY(N) f 8X»4HZ(N) »8X» 

14HU(N) ^8Xf 4HV(N) »8 Xf4HW(N) ) 

112 FORMAT{lX»I4rA(.lXrE11.4) ) 

113 F0RMAT{3X»6HTIME =»F9.3) - - . 

114 .FORMAT (/fZXf IH-I r 2Xt 4HC IRC f QXr 3HWvL f 7X 1 4HAMPX.f 7X 1 4HAMP Y f 
1 7Xf5HRCOREr5XflOHSPAN., L.OC. f IXr lOHDEL. . PHASE) 

115 FORMAT (30H NUMBER TOO LARGE FOR MAG TA^E) 

116 FORMAT (31H NUMBER TOO LARGE FOR CONT CALC) 

117 FORMAT (4F10.0) 

118 F0RMAT{/,»2Xr5HiPAiR»5Xf4HciRC) 

119 F0RMAT(2I4f4E13.6) 

120 FORMAT (3Fi0. Of 14) , - 

202 F0RMAT(2XflHNf5Xf2HXFf9Xf2HYFf9Xf2HZFf9Xf2HDXf 9Xf2HDVf9Xf2HD2f9Xf 
12HDR) 

203 FORMATClXf 7HNPAIR “fI4) 

300 F0RMATC3F10.4) 

301 FORMAT (IHl) 

302 FORMAT ( 20 A4) 

303 FORMAT (IHI) 


•44 



531 

532 

533 
53^1 

535 

536 


304 F0R^lAT(2Xf f4-HXCG=fE11.4f6H YCG=,E11.4f 9H VDRIFT=rEll*4f 

17H DYCG=fEll.i|) , . ■ 

305 FORMAT { 2X r 5HXMAC= r FI 0 . 4 ) 

998 CONTINUE 

STO° 

END 




1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22. 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 
34' 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 
-55 

56 

57 

58 


C 


PREDICTOR - CO^^RECTOR TIME- STEPPING ^ ^ • 

SUBROUTINE PCOR(Nf IPCO/IPr I8R» NPAiRr IHl r IHALFP? D£LT» iWr ISTf 10 r 

1 JAC t XMAC ) ~ ' , ^ .. . „ . . ‘ 

DIMENSION XBARC7) »YBAR(7) »ZBMI (7) rZBPL(7) ? ZBMI2,(7) ?:Z8PL2(7) 
COMMON X(92»7) rY(92f7) »Z(92f7) fU(92»7) »V(92r7) fW(92^7) f 

l'XDUM{93»7> f YDUM<93»7) f ZDUM(93f7) »CX{92r7) fCY(92r7) >CZ{?2f7) r 

2 CIRC (7) rALMBDA(7) » XO ( 92f 7) f YO ( 92f 7) » ZO C 92? 7) /DELRM{92»7) 

3 TRC0'M(92f7J 

DATA PI4/12, 566370/ 

Xlf^=I>3 


IF(IHALFP.LE.l) GO TO 12 
DO 5 U=lfNPAIR 
IF(IO.GE.i) GO TO 6 


X8AR{L) = {XDUM(2fU+XDUM{IHALF“»U )/2.0 
YBAR{L)={YDUM(2»L)+YDUM( IHALFPrL) )/2.0 
Z8M I C U =-XiP+ AUMBD A ( U /2 . 6 
ZB-^L ( L) = C xi'=-i . 0 ) yALMBDA ( L) /2. 0 

GO TO 5, . . 

6 XBAR(L) = tXDUM{2fU)+XDUM( IBR+3,fUl )/?.0 
YB AR ( U =.( YDUM ( 2 r L) + YDUM ( I BR+1 » L) ) >2 . 0 
IF,( JAC.LE.p) GO ,TO 13 
XB AR { U ) =xbUM ( IH ALFP t L ) 

YBAR { L) =YDUM (.IHALFPtL) 

13 ZBMI { U =- ( XIP-l. b ) YALMBDA { U 
ZB°L{L)={XIP-i.5) YALMBDA(U) . 

ZBMI2 C U =“ ( Xio-1 , 5 } Y ALMBDA ( L ) 

ZBpL 2 ( U = ( X I°-l . 0 ) Y ALMBDA ( L ) 
iF(JAC.LE.b) GO TO 5 
ZBMI2 ( U'=Z ( 2 f U - 6 . 375YXMAC 
ZBPL2(D=ZDUM( IHALFPfU " 

5 CONTINUE 
12 CONTINUE 

□SL^Oji^ 

— — VSL=b. 

wsl=6. 


DO 11 J=1»NPAIR 
DO iO LU=2fIHi 


L=LL 

IFdHALRP.LE.D - L=LLrl . 

CALL FAAyEL(N»i» Lf IPfUL> VLf WLf IBR» IHALF‘=^ IW» IST» 10) 
OSL=USL+UL 


VSL=VSL+VL 

10 WSL=WSL+WL , , , ' 

IF( IHALFP.LE.*,!) ,G0 TO 11 
iF(IO.LE.b) ,Gp TO i 
IR(JAC.GE.l) GO TO 14 
IF(N.UE.IBR) GO TO 1 
14 ZPL=ZBPL2(J) 

ZMi=ZBMI2(J) 

GO TO 2 

1 ZPL=ZBPL(J) 

ZMi=ZBMnU) 

2 DXB=XBAR( J)-X{NFi) 
DYB=YBAR'( J)-Y(Nd) - 
HSQ=0 XByDXB+DY8ydVb 
H=SQRT(HSQ) 

DXBi=-XB-AR ( O ) -X { N r I ) 
DYB1=DY8 
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59 


HSQl=DXBl*DXBl+DYBl’t=DYBl 

60 


H1=SQRT(HSQ1) 

61 


DZ^=ZPL-Z(Nfi) 

62 


DZM=ZMI-Z(Nf I) 

63 


IF(HSQ-.OOOi)3r4^4 

64 

3 

vel=6.6 

65 


GO TO 7 _ , 

66 

4 

COSP=DZP/SQRT ( DZ°^bZP4-HSQ) 

67 


COSM=DZM/SQRT ( OZM+DZC-l+RSQ) 

68 


IF{ JAC..GE.1) C0SM=”1^6 

69 


VEL=-CIRC( J)’f={2.0+(COSM-COSP) )/(PI4+H) 

70 

7 

lF(HSQi-.0D0i)8»9f9 

71 

8 

VEL1=0 , 0 

72 


GO TO 15 

73 

9 

C0SP1=DZP/SQRX ( DZ^*DZ‘='+HSQi ) 

74 


COSMl=OZM/SQRT ( DZM+DZM+HSQl ) 

75 


iF(JAC.GE.l) COSM1=-1,0 

76 


VELi=CIRC( J)>!=(2.6+(COSMi-COSPl) )/{PI4’i=Hi) 

77 

15 

lF(HSO“.6001)16fl7f 17 

78 

16 

C05A=i.O 

79 

* 

SiNA=1.0 

80 


GO TO 18 

81 

17 

COSA=DXB/H 

82 


S1NA=DYB/H 

83 

18 

IF(HSQi“.0001)19f26>26 

84 

19 

C05A1=1.0 

85 


SINAi=1.0 

86 


GO TO 21 

87 

20 

C0SA1=DXB1/Hi 

88 


SINA1=DYB1/H1 

89 

21 

USL=USL-VEL+SINA--VELl>»=siNAi 

90 


VSL=VSL+VEt*C0SA+VELl^C05Al 

91 

ii 

CONTINUE 

92 


U(Nf I)3:USL 

93 


V(NrI)=VSL 

94 


W(Nf i)=WSL 

95 


RETURN 

96 


END 
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C 


C 

C 

C 


t 


INDUCED VELOCITY SUBROUTINE ^ 

SUBROUTINE "FAAVEL (Nf II f 12^ I f'MMAXf US» VSf WS? IBRr IHALFPf IW? IST» I0P2) 
COMMON X(92t7) f Y(92f7) »Z(92»7) »U(92f7) fV{92f7) f W(92f7) ? 
i"XbUM(93»7) rYDUM(93f7) f ZDUM(93r7) rDCX(92>7) »DCY(92r7) ?DCZ(92r7 ) t 

2 CIRC(7) ?ALMBDA(7) ^X0(92f7) f Y0(92r7) ►Z0(92r7) j'DELRM { 92» 7) 

3 fRC0Mi92f7) 

data PI/3.1415926/fPI4/12. 565370/ 

ZM1=Z( If 12) 

Z|^1=Z,(N/I1) 

SIGN2=1.5 

iab?6 

ICD=0 


10 

20 

1 

6 

7 


8 

40 


36 


50 


US = 0;.0 
VS=0.0 
WS=6.0 


M=lfMMAX LOOP PICKS ^UP SUCCESSIVE^ HALF-WAVES.. 
EFFECTS OF BASIC VORTICIES PLUS REFLECTIONS IN THE 

other three quadrants 

CX=DCX( Jf 12) 

Cy=DCY(If 12) 


CZ=DCZ(If 12) 


XRATI=X(If 12)-X(Nf ID 
YRATl = Y { i 1 12 ) - Y ( N f ri ), 
XRAT2=-X(If I2)-X(Nf II) 
YRAT2=- Y U f 12 ) - Y C N f I i ) 
XRlS^XRATli^XRATi 
YRlS=YRATl*YRAtl 
XR2S=XRAT2^XRAT2 
YR2S=YRAT2*YRAT2 
DO 990 M=lfMMAX 


IM=MMAX-M 


U£=0.0 

V£=0.0 

WE-O-.O 


IF(MMAX,£Q.l)GO- TO 50 
IF( (l-M)+iM)^40fidf20 
IF{IM) 50flf50 
WRITE(IWf200) 

GO TO 999 
iF(I0P2)6f6f7 
ICD=i 
GO TO 40 


IF(I)i.GT.IBR) GO TO 8 
icD=i 


IFCI.LE.IBR) GO TO 32 

GO JO 40 

iAB=l 

IFd.GT^IBR) ,GO TO 32 
IF(iOP2.G£.l)GO to 36 
CX=-CX 

cY=-cy .. 

SIGN2=-SIGN2 

SIGNl=-(M/2) . ... 

Z(If i2)=SIGNl*ALM8 DA(J 2)+SIGN2+=ZMl _ 

IF ( 105=2. G£. 1 ) Z ( I f 12 ) =- < M-1 ) *ALM8DA ( 12 ) +ZM1 

ZRATl=Z(if-I2J-ZNl 

IF ( IHALFP. L£. 1 ) ZRAT1=0 . 

ZR15=ZRAT1+ZR,AT1 
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60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 
117 


IF(IAB.GE.i) GO TO 2 
RSQA=XRlS+YRlS+ZRiS 
IF(iOP2.LE.O) GO TO 37 , 

2 ZRAt2=Z ( I f 12) + ( 2*M-2) ?ALM8DA ( 12 ) ~ZN1 
IFCIHALFP.UE.l) ZRAT2=0. 

GO TO 38 

37 ZRAT2=-Z(I^I2)-ZN1 

38 ZR2S=ZRAT2’f=ZRAT2 
!f(Iab.g£»i) go to 3 
RSQB=XR2S+YRlS+ZRiS 
IF{IHAUF°,L£.l) GO TO 43 
IF(ICD.GE.l) GO TO 5 

3 RSQC=XR2S4-YR1S+ZR2S 
RSQD=XR1S+YR1S+ZR2S 

43 IF{IA8.G£.l) GO TO 4.. 

5 IF{RSQA-.0001)ll»12fl2 

11 R£CIPA=0.0 

GO 'JO 13 ,, 

12 telPA=l,.0/{pI4>t^RSQA?SQRt(RSQAn , 

i F C IHALFP . L£ . 1 ) RECIPA=1 . 0 / ( 2 . 0 *P I *RSQA ) 
i‘3 iF(RSQB“,000i)14»15f 15 

14 RECiP8=6,6 

GO “TO 16 , . . 

15 RECIPB=1.6/(PX4*RSOB%=SQRT(RSQ8n _ 
IF(IHAtFP,LE.i) RECIPB=1.0/(2*6+°I>t'RSQB) 

16 IFCICD.GE.i) GO ,TO '22 _ 

IF{ IHALFP. LE.i) GO TO 17 

4 IFCRSQC-,0001)17ri8fl8 

17 R£CiPC=0.b 

G0~t0'19. , . , 

18 RECI°C=1.0/(PI4*RSQC*SQRt(RSQC) ) 

19 IF(IHALFp.LE.i) GO TO 23 ' ‘ 

iF(RSQD-.000i)23f21f21 

23 R£CI‘=D=0.0 
GO 'TO 22 

21 R£CIPD=1.0/(PI4tRSQD+SQRT{RSQD) ) 

22 iF(IAB.,GE.l) GO TO .30 
UE=R£CIPA* { CZ+YRAT1-CY=>:ZRAT1 ) 

VE=RECIPA* ( CXi=ZRATlrCZ*XRATl ) 

W£:=R£CIPA’f= ( CY+XRATl-CX+YRATl ).. . 
U£=U£+R£CI°B’t‘ ( CY+ZRATl-CZ+YRAXl ) 
V£=V£+RECIpB^= { CX*ZRATl+CZ’»'XRAt2) 
WE=WE+RECIPB* ( -CY*XRAT2-CX+^YRAT1 ) 
IFUCD.GE.i) GO'TO 32 
iF(I0P2.GE.l) GO TO 30 . , 
U£=U£+RECIPC*{-CZ+YRAT1-CY*ZRAT2) 
VE=V£+RECiPC+(-CX*ZRAT2+CZ*XRAT2) 
V/£=WE+R£CipC* ( CY+XRAT2+CX+YRAT1 ) 
UE=UE+RECIPD4= { CZ*YRATH-CY*ZRAT2) 
V£=V£+RECIPD* ( “CX+ZRAT2-CZ+XRAT1 ) 
WE=W£+R£Ci°D* ( -CY+XRAT1+CX*YRAT1 ) 

GO TO 32 

30 iF(M.EQ.i) GO TO 32 

U£=U£+R£CiPC* { -CZ*YRATH-CY^=ZRAT2 ) 
VE=V£+RECIPC+ { +CX’^<ZRAT2+CZ*XRAT2) 
W£=WE+RECIPC+(-CY4=XRAT2“CX*YRAT1) 
UE=U£+R£C i PD 1= ( CZ^ YRATl-C.Y* 2RAJ2 ) 
VE=VE+R£CipD=f< { +CX+ZRAT2-CZ+XRAX1 ) 
W£=WE+R£CIPD’t= ( +CY+XRAT1-CX*YRAT1 ) 
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iia 

119 

120 
121 
122 
123 
124’ 

125 

126 


32 US=US-UE 
. VS=VS-V£ 

990 WS=WS~W£ 

Z(Ir I2)=ZM1 
. Z(Nril)=ZNl 

‘ 999 RETURN ^ ^ 

200 F0RM/\T(5(/) »7iH+*+DIAGN0STIC+**— — THE NUMBER OF PERIODS (IP) CANN 
SOT BE GREATER THAN 5) “ 

END ' ' 
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24 
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28 

29 

30 

31 

32 

33 

34 

35 

36 
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C CURVATURE EFFECTS SUBROUTINE 

SUBROUTINE CURV ( N t DELT^ IH ALF f IPCO 1 1 f f RNC f I p UC f VC f WC f I W f 1ST ) 

.COMMON, X( 92c' 7) cY{92c7) cZ(92c7) cU(92c75 cV{92c7) cW(92c7) c 
l'XDUM{95c7) f YDUM(93c7) c ZDUM ( 93c 7) c CX ( 92c 7) c CY { 92c 7) c CZ( 92c 7 ) c 

2 CIRCC7) cALMBDA(7) c x6 { 92c 7) c YO ( 92c 7) cZ0(92c7 ) cD£LRM(92c7) 

3 TRC0M(92c7) 

DATA 1415926/ 

itHET=0 

C RADIUS OF CURVATURE CALCULATION 

.ASQ={X(N+lci)-X(N-;ic“I) ) *<=2+ ( Y( N+1 c I ) -Y ( N->1 c I ) )*^2+ 

1 (Z[N+ic,,I)-Z(N~ic I) )*>K2 
ASSORT (ASQ)., „ ^ 

8SQ=(X(N+ici)~XCNcI) )*t= 2+( Y(N+lc I)-Y(Nci) ) ^*24- ( Z( N-fl c I ) -Z { N c I ) ) +*2 

B=SQRT(BSQ) . , . . 

CSO=(X(Nc I)-X(N-ic i) } + +^2+< Y(Nc I )-YCN-lc I) ) *i^2+(2(Nc I)-Z(N-lc 1 ) ) r^'2 
C=SQRT(CSQ) 

DELRSB+C 

CCX=( Y(Nc h-Y{N-lc I) U{Z(N4lcI)-Z(Nc I) )- ( Z( Nc I ) -Z( N-l c I ) ) * 

KYCN+ic I)-Y(NcI) ). 

CCY=CZ(NcI)>-Z{N-lci) )>«=(X<N+icI)“X{NcI) )“{X(NcI)-X(N-lci) ) + 
l(Z{N+i cl)::Z(!^c I) ) , 

CCZ= ( X ( N c n -X ( N-1 c I ) ) + { Y ( N+1 cI)-YCNcI))-(Y(NcI)~Y( N-1 c I ) ) t 
KXCN-J-icD-XCNcI) ) 

CCSQ=CCX!i=CCX+CCY*CCY+CC2+'CCZ 
iF(CCSO)4lc'5c5" " 

5 IF( (CCSQ/[CSQ>s=BSO) )”.00QCl)6c6c7 

6 iF(ASQr(BSQ+CSQ) )32c32c33 
32 COSA=6*707 

COSBs 6.7Q7 
COSCsO.O 
GO TO 22 
35 COSA=i,0 
COSBsi.O 
C0SC=1,0 
GO TO 23 

7 CC4sQRT(CCSQ) 

COSAsCCX/CC 

COSBSCCY/CC 

COSCrCCZ/CC 

iF(C)l3c42cl3 

13 i?:(B)llf42cll 
11 SINAL=CC/(C+B) 

IF ( ASQ- ( BSQ^CSQ ) ) 3 c 4 c 4 

3 QALRHA=1.0 ' 

GO TO 8 

4 QALPHA=2.0 

8 lF(SiNAL-,0oi)9c2c2 

9 rF(QAL°HA-l*5)22c22c23 
2 RC=A/(2.6*SINAL) 

GO, TO 26 

22 RC=(B+C)/2.0 
THl4=PI/2.b 
TH24=TH14 
ITHETsl 

GO TO. 12 

23 VELC=0.0 

RC0=SQRT{ RCOM{ Nc I ) *4=2+4. 0+DELT/RNC) 

ITH£T=2 
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GO TO 12 , 

C :CORE RADIUS AND SELF INDUCED VELOCITY CALCULATION 

' 26 IF(B-2.0’i‘RC)28r43r43 ~ 

28 iF{C72«0+RC;27»43r43 

27 THi4=6.25+ATAN(B/(SQRX(4.6+RC’t=RC”B+8) ) ) 
TH24=b.254=AtAN(C/(SQRT(4.6*RC^'RC-C=HC) ) ) 

12 IF(IT.GT»i;GO TO 34 - ~ 

IF ( I PCO . EQ . 6 ) DELR.M ( N r I ) =DELR 
34 IFaTH£T,E'a,2) GO JO 10 
16 D2=RC0M(I::JrI)>t=+2*{DELR-DELR.M(Nr i) ) /DELR 
DRC02=4. 6+DELT/RNC-D2 
RC02=RC0M(NrI) *=»=2+DRC02 
, iF(RC02)44r44f29 

29 RC0=SQRT(RC02) 

. IE(RC)45r45r30 

30 DV1=ALOG(8.0+RC/RCO)~0.558 
IF(ITHET)18fl8fl9' 

19 qV2=0.0 

GO to, 20 , , 

18 DV2=0.5*ALQG( (C0S(THi4)/SJN{JH14.) )^{C0S(TH24)/sfN(TH24) ) ) 

20 V£LC=CIRC(I),*(bVl-DV2)/{4.6+°I’t=RC) " 

10 UCtrVELC^^COSA 

VC=VELC*COSB 
WC=V£LC*COSC 
rF(iPCO.GE.I,)GO to 1 
D£LRM(Nti3=D£LR 
RCOM(Nf I)=RCO 
GO^TO i 

41 WRIT£(iW»ibO) 

GO, TO 46 

42 WRltE(IWflOl) 

GO TO 46 

43 WRitE(iWfid2) 

GO, TO 46 

_ „44 

GO to 46 

45 WRIT£(iWfl04) 

46 ISt=l 

1 CONTINUE 

RETURN , . 

100 FORMAT (15H AB5,(C>;8) .UT,0* ), , 

101 F0RMd(23H AQJ» POINTS COINCIDENt) 

102 FORMAT (24H RAD, OF CURV, TOO SMALL) 

103 F0RMAT(i2H RCORE.LEIO. ) , . , 

104 F0RMAT(19H RAD. 'OF‘CURV.Le.6, ) 

END 
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C COMPARISON Of GEOMETRY- AfTER BK'/iD. MARCH WITH INIT, GEOMETRY 

SUBROUTINE STAT ! Nf I r‘ IHI » XAMf XGMe IWJ T . . ‘ 

COMMON X(92r7) ( Y(92>7) ^ Z(92f 7) rU{92ii7) f V(92f 7) e W(92»7) f 
1~XDUM(93>7) fYDUM{93f7) f ZDUM{93f 7) »CX(92f 7) r CY (92 f 7) f CZ{ 92 r 7 ) t 

2 CIRC(7) f ALMBDA(7) f x6 ( 92^ 7) f YO ( 92t 7 ) r Z0(92>7) >D£LRM(92f 7) 

3 rRC0M(92r7) 

DX=X(N»I)-XO(Nf I) 

DY=Y(Nf I)-yO(Nf I) 

DZ=Z<N^I)-z6(N»I) 

DR=SQRT ( DX+DX+-DY*DY-fDZ^=DZ ) 

XNP=N 

XAM=XAM+DR 

XGM=XGM*DR 

NP=N-1 

WRiT£(6^200)NPf XCNf I) > Y ( N f I ) » Z( N f I ) f DX» DYf DZ » DR 
If(N.LT.IHl) GO TO 25 
XAME=XAM/XNP . 

XGME=(XGM)*^'(1./XNR) 
write ( 6^201) XAME r XGME 
25 CONTINUE 
RETURN 

200 FORMAT(i3f7Eli.4) .... 

201 F0RMAt(2X»16HARiTH MEAN DEV E11.4r2Xf iSHGEOM .MEAN DEV =fE11.4) 

END 
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Vn 

■vj 


NOTE; DISTANCES NORMALIZED ON WING SEMI -SPAN 
Figure 2.1 Two-Dimensional Twin Vorte>< Pair Trajectories 


VERTICAL VEL. /DRIFT VEL 







INITIAL CONDITIONS 

r = 0.127 u b 

CO 

X = 5.0 b 

= 3y = 0.125 b 



CALCULATION FROM REF. 
10 POINTS/HALF WAVE 
CORE UNSPECIFIED 
TIME = 5.33 b^/r 



PRESENT CALCULATION 
10 POINTS /HALF WAVE 
CORE RADIUS = .098 b 
TIME 5-33 b2/r 


Previous and Present Calculations of the Distortion of 
a Vortex Pa i r 



.TIME-TO-LOOP PARAMETER rt( 



PERTURBATION AMPLITUDE, ax/b, Sy/b 


Figure 3-2 Dependence of Time-to-Loop upon Initial Perturbation 
Amplitude, X= 5.0, 6/b = .098 
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TIME-TO-LOOP PARAMETER 'Tto/h^ 





8 ^ 


4 -^ 





0 L 
0 


T 

.1 


VORTEX CORE RADIUS 5/b 


.4 


Figure 3,4 Dependence of Time- to- Loop Upon Vortex Core Rad 
a^/b = a„/b = 0 , 125 , X = 5.0 



TIME-TO-LOOP PARAMETER Ftj/b^ 



Figure 3-5 


Dependence of Time-to-Loop Upon Perturbation Wavelength, 
a^/b = ay/b = .125, 6/b = .098 
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VERTICAL DISTANCE/b 



Figure 3.6 End-View of Vortex-Pair Distortion Toward Looping, 

X = 5.0, d/b = 0.098 


6k 


VERTICAL DiSTANCE/b 





V|_H -STRAIN VELOCITIES DUE 
TO L.H. VORTEX (DRIFT 
REMOVED) 

~ SELF- INDUCED VELOCITY 
DUE TO CURVATURE 

Vres ~ resultant of above 


V^^f^p -- VELOCITY TENDING TO AMPLIFY WAVE 

. VELOCITY TENDING TO ROTATE WAVE 'BACKWARDS' 


Figure 3*8(b) Mechanism for Oppos i ng-Sense Vortex and Wave 
Rotations - Velocities 














VERTICAL DIST 







VERTICAL DIST 



SPANWISE DIST 
(a) ry/rp = 0.6 

NOTE: VERTICAL SCALE COMPRESSED 

BY A FACTOR OF 5- 



SPANWISE DIST 
(b) Fy/rp = 1.0 



SPANWISE DIST 
(c)^ Fy/Fp = 1.667 


Figure Cross-sectional Trace for the Convergence Point on Twin Vortex 

Systems as a Function of Fy/Fp (np = 0.333, Hy =1.0) 


<uO 




SPANWISE 

POSITION 



SPANWISE 

POSITION 



Figure h.5 Examples of Tip- and Flap-Vortex Convergence 
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O FLAP CONVERC 
A TIP CONVERGE 
EQUAL CHANCE 


Figure h.S Twin-Pair Parametric Study: 
Convergent Vortex 


Identification of the 
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Figure ^.9 Strength of the Vortex Remaining After Convergence 






SUCCESSIVE 'FRAMES' 
OF THE CALCULATION 


Figure ^.12 The Application of a Space-Fixed Calculation 







DIST. BETWEEN VORTICES AT t = 50 (2-D CALC.) 

DIST. BETWEEN VORTICES AT T = 0 

EFFECTS OF REFLECTION ACROSS CENTER PLANE 
NOT INCLUDED 



<3> 


<$> A 


. 


STEPPING INDEX.(A0) 




Figure Effect of Integration Step-Size on Vortex Spacing 
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30(3) 


— 

SAWTOOTH 

1 1 1 1 1 1 * 


0 O2O<iO6o 8I0 
(a) CV990, SAUTOOTM 


0.3H 
0 2 
0 H 
0 


23(2) 


10, 15(3) 
lOG 


0 02 0 l| 06 o 8 1.0 

(b) B7'i7, FULL CONFIGURATION 
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HUNS 2, 3. 1), 5, 29 2-D INITIALIZATION . RUNS 9, 10, II, 12, 2*1, 3*0 RELAXATION, A/C FIXED 

RUNS I'l. 15. 16. 17 3 FAIRS, SPACE FIXED RUIIS 22. 23. 2l|. 25 2 PAIRS, SPACE FIXED 

FIGURES IN PARENTHESES SHOW NUMBER OF PAIRS 
FIGURES WITHOUT PARENTHESES SHOW RUN NUMBER 


Figure 5.1 Span Load and Run Definition for Multiple-Pair Cases 
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NEAR WAKE 
{SECTION 5) 

o MULTIPLE PAIRS 
o VORTEX MERGING 
o VORTICES CONVECTED AWAY 

COMPUTING 

o AIRCRAFT-FIXED RELAXATION 


MIDDLE WAKE 
(SECTION k) 


FAR WAKE 
(SECTION 3) 


o TWIN PAIR REMAINS 
o SPACING CHANGES DUE TO OTHER ' 
VORTICES 

o POSSIBLE DISTORTION OF SPIRAL, FORM 
COMPUTING 

o AIRCRAFT-FIXED, FROM NON-SPIRAL FORM 
o SPACE-FIXED FROM PARAMETRIC, 

IF NEAR-SPIRAL 


o SINGLE ORGANIZED PAIR MAY REMAIN 
o LOOPING POSSIBLE 

o DECAY DEPENDS ON PREVIOUS DISTORTION 
COMPUTING 

o SPACE-FIXED FOR "CROW" INSTABILITIES 
o VISCOUS ANALYSIS FOR BURST OR SLOW VISCOUS DECAY 



NOTE. IN THE ABOVE SKETCH THE AXIAL DIMENSION IS FORESHORTENED AND THE 
NUMBER OF WAVES IS REDUCED. 


Figure 5*2 Calculation Regions for a Mui tipie-Pai r Vortex Wake 
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(a) Near-Cyclic Configuration 


TIME = 40 
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(b) Trans'ient Configuration 


Figure Neai — Cyclic and Transient Cases 
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NOTE- "TICKS ON AXES REPRESENT ONE SEMISPAN VERTICALLY 
AND SPANWISE AND ABOUT SIX SEMISPANS AXIALLY 


Figure 5»7 Effect of the Outboard Flap and Landing Gear on the Near V^ake 
of a Boeing ']h^ on Approach 


